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Abstract 


Numerical  modeling  of  shock  propagation  and  reflection  is  of  interest 
to  the  Department  of  Defense  (DoD).  Propriety  state-of-the-art  codes  based 
upon  E.  F.  Toro’s  weighted  average  flux  (WAF)  method  are  being  used  to 
investigate  complex  shock  reflection  phenomena.  Here  we  develop,  test,  and 
validate  a  one-dimensional  hydrodynamic  shock  code.  We  apply  WAF  to 
Gudonov’s  first-order  upwind  method  to  achieve  second-order  accuracy. 
Oscillations,  t5^ical  of  second-order  methods,  are  then  removed  using 
adaptive  weight  limiter  functions  based  upon  total  variation  diminishing 
(TVD)  flux  limiters.  An  adaptive  Riemann  solver  routine  is  also  implemented 
to  improve  computational  efficiency.  This  one-dimensional  code  is  then 
extended  into  two  dimensions  via  Warming  and  Beam’s  variation  on 
dimensional  sphtting.  The  numerical  capabihties  of  the  two-dimensional 
code  are  demonstrated  by  modehng  the  detonation  of  a  cylindrically 
symmetric  explosive  with  the  axis  of  the  cylinder  oriented  horizontally  above 
an  ideal  surface. 
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Chapter  I:  Introduction 


A.  Motivation 

Numerical  modeling  of  shock  propagation  in  hydrodynamic  flows  is  a 
field  of  obvious  interest  to  the  Department  of  Defense  (DoD).  Some  aspects 
are  particularly  difficult  to  model,  for  example,  the  transition  from  regular  to 
Mach  reflection.  Recent  state-of-the-art  codes  based  upon  E.  F.  Toro’s 
weighted  average  flux  (WAF)  method  are  being  used  to  investigate  such 
shock  reflection  phenomena.  Unfortunately,  these  codes  are  company- 
proprietary  and  unavailable  for  use  within  DoD.  The  objective  of  this  effort  is 
to  develop  our  own  WAF  code  for  use  at  the  Air  Force  Institute  of  Technology 
and  throughout  DoD.  Ideally,  we  will  find  improvements  in  accuracy  and/or 
computational  efficiency  in  comparison  with  the  algorithms  used  in  other 
codes. 

B.  Background 

In  1959,  Godunov  introduced  an  extension  of  the  first-order  upwind 
Courant-Isaacson-Rees  scheme  for  solving  the  Euler  Equations  (LeVeque, 
1992).  This  scheme  was  the  first  conservative  upwind  technique  developed 
and  has  given  rise  to  many  of  the  modern  numerical  techniques  in  use  today. 
Godunov’s  method  assumes  that  a  piecewise  constant  distribution  of  the 
conserved  variables  across  the  computational  domain  may  be  treated  as  a 
series  of  local  Riemann  problems.  The  flux  is  calculated  by  solving  the 
Riemann  problem  exactly,  forward  in  time.  While  robust,  this  technique  is 
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computationally  inefficient  because  each  Riemann  problem  requires  iteration 
to  solve  nonlinear  equations.  In  an  effort  to  make  this  scheme  more  efficient, 
several  approximate  Riemann  solvers  have  been  developed  that  solve  the 
Riemann  problem  without  requiring  iterative  processes.  Modern  Godunov- 
type  methods  are  constructed  in  this  fashion. 

First-order  schemes  are  known  to  result  in  solutions  smeared  at 
discontinuities  due  to  numerical  viscosity  and  Godunov’s  method  is  no 
exception.  In  1989,  E.  F.  Toro  presented  a  new  high-resolution  technique 
that  fully  utilizes  the  wave  structure  of  the  Riemann  problem  to  calculate  a 
more  accurate  intercell  flux.  This  flux  is  a  weighted  average  of  the  flux 
vectors  across  the  Riemann  problem.  Second-order  accuracy  is  achieved  in 
this  way.  Spurious  oscillations  t5q)ical  of  second-order  schemes  are  removed 
using  total  variation  diminishing  (TVD)  weight  limiters  (Toro,  1992). 

Extension  into  two  dimensions  is  easily  achieved  via  dimensional 
splitting  (Warming  and  Beam,  1976).  In  this  approach,  the  two-dimensional 
initial  value  problem  (IVP)  solution  is  obtained  by  solving  a  sequence  of  one¬ 
dimensional  problems  in  each  coordinate  direction.  This  technique  is  known 
as  time-operator  splitting. 

C.  Problem 

The  primary  objective  of  this  research  is  to  develop  a  two-dimensional, 
computationally-efficient,  hydrodynamic  shock  code  based  upon  the  weighted 
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average  flux  (WAF)  method  (Toro,  1989)  for  the  investigation  of  air  blast 
phenomenology. 

D.  Scope 

The  scope  of  this  project  is  limited  to  the  modeling  of  air  blast  under 
ideal  conditions. 

E.  Assumptions  and  Limitations 

The  computational  model  presented  here  assumes  that  air  behaves  as 
an  ideal  gas  with  a  constant  ratio  of  specific  heats,  i.e.  y  =  1.4.  This 
assumption  limits  use  of  the  code  to  low  pressure  regimes  where  peak 
pressures  are  approximately  690  kPa  or  below. 

F.  Approach 

Development  of  the  two-dimensional  code  begins  first  with  the 
introduction  and  application  of  all  numerical  techniques,  except  dimensional 
splitting,  in  one  dimension.  Code  development  is  performed  in  a  logical  step- 
by-step  manner  beginning  with  development  of  E.F.  Toro’s  exact  Riemann 
solver  and  Godunov’s  first-order,  upwind  method.  To  improve  computational 
efficiency,  an  adaptive  Riemann  solver,  also  introduced  by  Toro,  replaces  the 
exact  Riemann  solver.  Greater  stability  and  accuracy  is  achieved  near  steep 
gradients  using  an  adaptive  Courant-Friedrichs-Lewy  coefficient. 

Once  the  underlying  Godunov  algorithm  is  complete,  Toro’s  WAF 
method  is  added  to  obtain  second-order  accuracy.  Monotonicity  is  established 
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by  applying  one  of  four  TVD  weight  limiters  subsequently  written  into  the 
WAF  code. 

To  illustrate  each  numerical  tool  and  its  impact  on  computational 
efficiency  and  accuracy,  a  classic  one -dimensional  shock  tube  problem  is 
evaluated  over  a  course  mesh  for  comparison  with  the  solution  obtained  via 
the  exact  Riemann  solver.  While  this  solution  is  not  strictly  exact,  because 
no  analytical  solution  exists,  it  is  suitable  for  use  as  a  qualitative  benchmark. 
Once  the  one-dimensional  code  has  been  tested  it  will  be  vahdated  against 
experimental  results  obtained  from  the  Army  Research  Lab’s  57  cm  shock 
tube. 

Having  shown  the  one-dimensional  implementation  to  be  sound,  the 
code  is  then  extended  into  two  dimensions  via  the  dimensional  sphtting 
scheme.  To  illustrate  the  capabilities  of  the  code,  a  two-dimensional  problem 
is  solved  that  is  representative  of  the  types  of  shock  problems  of  interest  to 
DoD. 


4 


Chapter  2:  Governing  Equations  and  Important 
Thermodynamic  Relations 
A.  Simplifying  Assumptions 

The  basic  equations  of  fluid  motion  are  derived  from  the  conservation 
of  mass,  momentum  and  energy.  Full  consideration  of  these  conservation 
laws  results  in  the  Navier-Stokes  equations  of  fluid  dynamics  (Anderson, 
1982).  These  partial  differential  equations  (PDFs)  model  flow  with  a  high 
degree  of  fidelity  but  are  complicated,  computationally  inefficient  and 
impractical  for  many  apphcations.  The  solution  procedure  can  be  simplified 
significantly  if  the  following  assumptions  are  made: 

-  viscosity  is  negligible, 

-  heat  conduction  is  negligible  and 

-  gravitational  force  upon  air  molecules  is  negligible. 

For  propagation  of  shocks  over  an  ideal  surface,  the  inertia,  pressures 
and  physical  domain  involved  are  large  in  comparison  to  the  thin  boundary 
layer  and  its  effects  along  the  smooth  flat  surface.  This  means  that  the 
pressure  variations  normal  to  the  surface  are  negligible  through  the 
boundary  layer  and  can  be  ignored.  It  is  assumed  then,  that  the  pressure 
distribution  at  the  surface  is  due  entirely  to  the  flow  occurring  above  the 
boundary  layer,  therefore,  viscosity  can  be  ignored. 
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Heat  conduction  between  the  surface  and  air  can  also  be  ignored.  As 
stated  above,  an  ideal  surface  is  flat.  Another  important  property  of  our 
hypothetical  surface  is  that  it  is  perfectly  reflective.  Since,  heat  conduction 
can  not  occur  across  a  perfectly  reflective  surface  it  can  be  ignored. 

Lastly,  viscosity  and  heat  conduction  within  the  air  itself  can  be 
ignored.  Assuming  the  atmosphere  is  homogeneous,  the  velocity  and 
temperature  gradients  within  the  flow  field  are  very  small.  As  a  result 
viscosity  and  heat  conduction  within  the  air  can  be  neglected.  Under  such 
conditions  the  air  is  said  to  be  adiabatic  and  inviscid. 

An  important  thermodynamic  property  arises  from  these  assumptions. 
If  viscosity  and  heat  conduction  are  ignored  then  the  flow  is  considered  to  be 
isentropic.  This  means  that  as  a  particle  travels  through  the  medium,  the 
specific  entropy  at  the  particle  remains  constant.  The  nature  of  the 
compressible  flow  and  its  interaction  with  the  surface  is  now  fully  defined.  It 
is  inviscid  and  adiabatic  and  therefore  isentropic. 

One  important  exception  occurs  at  the  shock  fi*ont  where  compressive 
forces  result  in  significant  velocity  and  temperature  gradients.  These 
gradients  provide  the  necessary  dissipative  mechanism  for  the  transfer  of 
energy  and  momentum  across  the  shock  (Courant  and  Friedrichs,  1948). 
Here,  viscosity  and  heat  conduction  effects  are  significant,  the  specific 
entropy  is  no  longer  constant  and  our  simplifying  assumptions  are  in 
jeopardy. 
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These  effects  are  dispersive  and  result  in  a  thin  shock  front  across 
which  the  solution  is  smooth.  However,  the  width  of  this  region  is  roughly 
the  mean  free  path  of  the  air  molecules  and  is  microscopic  when  we  again 
consider  the  length  of  the  physical  domain.  Given  the  steep  gradients  that 
occur  and  the  relative  width  of  the  shock  front,  we  can  model  the  shock  front 
mathematically  as  a  discontinuity  and  once  again  ignore  viscosity  and  heat 
conduction  as  a  good  approximation.  In  fact,  if  we  obtained  the  smooth 
solution  using  the  Navier-Stokes  equations  and  compared  it  with  the 
discontinuous  solution  obtained  via  our  assumptions,  it  would  be  difficult  to 
distinguish  one  solution  from  the  other  (LeVeque,  1992). 

It  is  important  to  note  that  shocks  are  irreversible  thermodynamic 
processes  that  result  in  an  increase  in  entropy,  i.e.  specific  entropy  is  not 
constant  across  this  region.  In  modeling  the  shock  front  as  a  discontinuity, 
viscosity  and  heat  conduction  effects  are  not  ignored  but  are  approximated  in 
such  a  way  that  their  terms  can  be  eliminated  from  the  governing  equations. 
In  this  region,  the  flow  is  not  isentropic  and  therefore,  isentropic  relations  do 
not  apply. 

B.  Governing  Equations 

Application  of  the  assumptions  made  in  Section  A  above  reduces  the 
complex  Navier-Stokes  equations  to  another  coupled  system  of  PDEs  that  are 
nonlinear  and  strictly  hyperbolic  (LeVeque,  1992).  This  set  of  equations  is 
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commonly  referred  to  as  the  Euler  equations  or  hyperbolic  conservation  laws 
of  gas  dynamics. 

In  general  vector  notation,  these  equations  are: 

Ut  +  F{U)^  +  G{U)y  +  H{U),  =  0 .  ( 1 ) 

Note  that  each  vector  above  is  a  function  of  both  space  and  time  where  x  ,y , 
and  z  are  the  Cartesian  coordinates  and  t  represents  time.  Throughout  this 
document,  subscripted  independent  variables  indicate  partial  derivatives. 

In  Equation  (  1 ), 


C7  = 


P 

p  u 
pv  , 
piv 
E 


(2) 


is  a  vector  of  the  conserved  quantities:  mass,  momentum  (x  component), 
momentum  (y  component),  momentum  (z  component)  and  energy.  More 
accurately,  each  vector  element  above  represents  a  density  function  for  the 
corresponding  conserved  quantity  such  that 

\U^ix,y,z,t)dV 


is  the  total  quantity  of  the  m  conserved  variable  over  the  control  volume  at 
time  t  (LeVeque,1992).  In  Equation  (  2 )  above, p  is  density,  and  u  ,v ,  and 
w  are  the  components  of  particle  velocity  in  the  x ,  y  and  z  directions 
respectively.  E  is  the  total  energy  density. 
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(4) 


E  =  \pVV  +  pe, 

where  e  is  the  specific  internal  energy  to  be  defined  later  and  V  is  the 
particle  velocity  vector. 

Vectors 


F(U)  = 

p  u 

pu^  +  p 
pvu 

,  G(U)  = 

pv 

puv 

pv'^+p 

and  H(U)  = 

pw 

puw 

pvw 

pwu 

u(E  ^p) 

pwv 

v(E  +  p) 

pw^  +  p 
w(E  p) 

represent  the  flux  of  these  conserved  variables  in  the  x ,  y  and  z  directions 
respectively. 

Variables p,u  ,v,w ,  and p  are  often  called  primitive  to  distinguish 

them  from  the  conserved  quantities  (Toro,  1997).  Primitive  variables  are 
those  variables  that  can  be  controlled  and  measured  experimentally.  In 
practice,  it  is  convenient  to  transform  Equation  ( 1 )  into  its  primitive  form: 

W,  +  F(W)^  +  G(IV),,  +  H(W),  =  0 .  ( 6 ) 

where  W  can  be  expressed  in  terms  of  U  as 


u 

V 

= 

U3M 

w 

U4/C/1 

kP) 

,(r  -  i)(t/6  -  (0.5/  c;,  ){u2‘  \ 
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Equation  ( 1 )  can  be  simplified  further  if  we  take  advantage  of 
problem  symmetry  and  assume  uniform  flows  accordingly.  If  flow  in  the  0 
direction  is  uniform,  Equation  ( 1 )  reduces  to  the  two-dimensional  Euler 
equations, 

Ut+F(U)x+GiU)y=0.  (8) 

An  infinite  high-pressure  cylinder  is  an  example  of  this  type  and  will  be 
explored  in  Chapter  6.  Similarly,  if  we  assume  uniform  flow  conditions  in 
both  y  and  2 ,  reduction  of  Equation  ( 1 )  results  in  the  one-dimensional 
Euler  equations, 

Ut+F(U)^=0.  (9) 

Modeling  constant  cross  sectional  area  shock  tubes  is  a  classic  application  of 
Equation  ( 9 ).  In  Chapters  3  and  4,  problems  of  this  type  are  explored.  In 
either  case,  we  are  faced  with  a  dilemma.  There  are  m  + 1  unknowns  in  m 
equations,  so  the  problem  is  underspecified. 

C.  Important  Thermodynamic  Relations 

To  solve  the  Euler  equations,  an  appropriate  equation  of  state  is  used 
to  derive  a  closure  condition.  Here  we  assume  the  atmosphere  to  be 
comprised  of  a  diatomic  polytropic  gas  where  pressure  is  related  to  density 
and  temperature  by  the  ideal  gas  equation  of  state: 

p  =  pRT  ( 20 ) 
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where  R  is  the  specific  gas  constant  and  T  is  the  temperature.  For  an  ideal 
gas,  R  is  related  to  the  specific  heats  at  constant  pressure  and  volume  such 
that 

^  =  (11) 

where  Cp  is  the  specific  heat  at  constant  pressure  and  is  the  specific  heat 

at  constant  volume.  By  polytropic  we  mean  that  the  specific  internal  energy 
is  a  function  of  temperature  alone  such  that 

e  =  c^T.  (12) 

Direct  substitution  of  Equations  ( 11 )  and  ( 12 )  into  Equation  ( 10 )  yields 
the  following  important  thermodynamic  relation  for  specific  internal  energy, 
pressure  and  density: 


e  =  - 


(13) 


(r-i)p 

where  y  is  defined  to  be  the  ratio  of  Cp  to  c^.  For  air  at  pressures  below  100 


psi,y  wl.4. 


As  stated  earlier  in  Section  A,  the  flow  is  isentropic  everywhere  except 
at  the  shock  front.  This  leads  to  the  important  isentropic  relation 

pccp^.  ( 14 ) 

When  this  relation  is  applied  in  Equation  ( 10 ),  the  following  important 
isentropic  relation  is  obtained: 
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£_ 

Po 


\ 

P_ 


(16) 


where  the  subscript  zero  indicates  initial  conditions. 

Another  important  quantity  in  the  study  of  compressible  flows,  and  in 
particular  flows  where  shocks  occur,  is  the  local  speed  of  sound,  c ,  where 

^  dp\ 


(16) 


Differentiation  of  Equation  ( 14 )  yields  the  local  speed  of  sound  for  a 
polytropic  gas  given  constant  entropy: 


c  = 


(17) 


At  this  point  we  now  have  all  the  principle  equations  necessary  to  build  our 
numerical  methods  and  solve  shock  problems.  In  the  next  chapter,  we 
introduce  and  solve  an  initial  value  problem  (IVP)  of  importance  called  the 
Riemann  problem.  As  we  shall  see,  this  IVP  forms  the  backbone  of  all  our 
numerical  techniques. 
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Chapter  3:  The  Riemann  Problem  and  Its  Solution 

A.  The  Riemann  Problem 

The  Riemann  problem  for  the  one-dimensional  Euler  equations  is  the 
IVP  with  initial  conditions  (Chang  and  Hsiao,  1989): 


tj  (x,0)  = 


Ui(x,0) 

< 

U^(x,0) 


x<0, 

X  >  0. 


(18) 


The  domain  is  centered  about  the  discontinuity  at  x  =  Om  and  includes  all 
points  {x,t)  in  the  x-t plane  such  that  - oo  <  x  <  oo  and  t>0s.  In  practice, 
the  domain  is  restricted  about  a  finite  spatial  interval,  -  Ax/2  <  x  <  Ax/2 . 

An  important  mathematical  property  of  the  solution  is  its  self¬ 
similarity  nature  (LeVeque,  1992;  AGARD,  1961).  Provided  the  initial 
discontinuity  is  located  at  x  =  0/n ,  the  solution  is  a  function  of  one  variable, 
4'  =  xlt,  such  that 


Uix,t)  =  U(C).  (19) 

Given  Equation  ( 19 ),  the  solution  is  said  to  be  similar  at  different  times 
along  characteristic  curves,  where  4"  =  constant,  in  the  x  -  <  plane.  It  follows 

that  because  ^  =  constant,  dxj dt  =■  constant  so  the  solution  propagates  along 
these  characteristic  curves  at  a  constant  speed  as  well. 

This  self-similarity  property  forms  the  basis  of  a  powerful  technique 
used  to  solve  linear  first-order  hyperbolic  PDEs  called  the  method  of 
characteristics  (AGARD,  1961;  Abbott,  1966;  Anderson,  1982).  In  the  method 
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of  characteristics,  the  solution  to  a  linear  hyperbolic  system  of  m  PDEs  is 
obtained  by  solving  two  sets  of  m  ordinary  differential  equations  (ODEs):  set 
one  defines  the  characteristic  curves  along  which  the  solution  remains 
constant  and  propagates  and  set  two  defines  compatibility  equations  which 
hold  true  along  these  characteristics.  Unfortunately,  the  Euler  equations  are 
nonlinear  resulting  in  discontinuities  where  characteristics  cross,  so  the 
solution  procedure  is  not  so  direct.  Despite  this,  much  of  the  characteristic 
information  obtained  via  the  method  of  characteristics  still  remains  useful. 

B.  General  Wave  Structure  of  the  Riemann  Problem 

The  method  of  characteristics  plays  an  important  role  in  defining  the 
general  wave  structure  of  the  Riemann  problem.  Using  the  method  of 

characteristics,  we  obtain  three  characteristic  fields,  C~ ,  C®,  and  C^with 
characteristic  speeds  (AGARD,  1961;  Toro,  1997): 


=  u-c, 

(20) 

^2  =U 

(21) 

and 

Ag  =  U  +  C 

(22) 

respectively.  The  superscripts  -,  0,  and  +  reflect  the  role  of  c ,  the  local  speed 
of  sound,  in  the  above  equations. 

Each  characteristic  field  presented  above  defines  a  wave  in  the  RP 
domain  such  that  wave  1,  is  defined  by  C~  characteristics,  wave  2  by 
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C®  characteristics  and  wave  3  by  characteristics.  In  the  :r  -t  plane,  these 
waves  partition  the  RP  domain  into  four  regions  as  shown  in  Figure  1: 
Ri,R2,R^  andi?4  where 


U (x,  t)  has  been  recast  into  its  primitive  variable  form  using  Equation  (  7 )  for 

convenience.  Throughout  most  of  our  discussion  and  computations,  the 
primitive  variable  form  can  and  will  be  used. 

The  outer-waves,  1  and  3,  can  be  either  rarefactions  or  shocks 
depending  upon  their  characteristic  fields.  If  the  characteristic  field  is 
divergent  the  wave  is  a  rarefaction;  if  it  is  convergent  the  wave  is  a  shock  (see 
Figures  2  and  3).  Until  the  true  nature  of  the  wave  is  determined,  each  wave 
is  represented  by  a  pair  of  rays  meant  to  represent  the  head  and  tail 
characteristics  of  a  rarefaction.  The  inner-wave,  wave  2,  is  always  a  contact 
surface  and  is  represented  by  one  ray.  The  large  arrows  in  Figure  1  show  the 
motion  of  these  waves. 
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Figure  1:  General  Wave  Structure  of  the  Riemann  Problem 
If  the  outer  wave  is  a  rarefaction,  the  solution  varies  continuously 
along  divergent  characteristics  that  lie  between  the  head  and  tail.  For 

example,  we  see  that  the  C®  characteristic  field  defines  a  left  traveling 
rarefaction  wave  as  depicted  in  Figure  2.  Because  the  characteristic  field  is 

t 


head 


Figure  2:  Charcteristics  Defining  a  Left  Traveling  Rarefaction 
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divergent,  pressure  and  density  decrease  as  we  travel  across  the  expansion 
fan  from  head  to  tail.  The  particle  velocity  may  increase  or  decrease 
depending  upon  whether  the  rarefaction  is  traveling  to  the  left  or  right.  If 
wave  1  is  a  rarefaction,  particle  velocity  increases.  If  wave  3  is  a  rarefaction, 
particle  velocity  decreases.  The  speed  of  the  rarefaction  wave  is  taken  to  be 
that  of  the  head  characteristic. 

If  the  outer  wave  is  a  shock,  characteristics  converge  and  cross  as 
shown  in  Figure  3  for  a  right  traveling  shock.  The  solution  is  discontinuous 


t 


Figure  3:  Characteristics  Defining  a  Right  Traveling  Shock 
at  the  point  where  they  meet.  By  convention,  the  head  and  tail  are  dropped 
and  the  shock  front  is  represented  by  a  bold  solid  line  that  represents  points 
where  these  characteristics  cross.  Each  primitive  variable  experiences  an 
instantaneous  jump  across  this  line  as  we  cross  the  pair  of  rays  from  head  to 
tail.  The  speed  of  the  shock  wave  is  not  that  of  the  head  characteristic  since 
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this  wave  itself  is  not  a  characteristic.  The  speed  of  the  shock  wave  will  be 
introduced  later. 

As  stated  above,  the  middle  wave  is  always  a  contact  surface.  Here, 
characteristics  run  parallel  to  the  contact  surface  as  shown  for  a  right 
traveling  contact  surface  in  Figure  4.  Particle  velocity  and  pressure  are 
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Figure  4:  Characteristics  Defining  a  Right  Traveling  Contact  Surface 

constant  across  the  contact  surface,  i.e. 

W.  =  M2=«3  (24) 


and 


P*=P2=Pi-  (25) 

Because  there  is  no  pressure  differential  across  wave  2  and  particle  velocities 
are  the  same  on  either  side,  the  flow  remains  separated  resulting  in  a 
discontinuity  in  density.  The  region  that  lies  between  waves  1  and  3  is  called 
the  star  region,  .  Note  that  f?*  consists  of  i?2  -^3  • 
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The  complicated  nature  of  the  waves  described  above  means  the 
solution  of  the  Riemann  problem  can  have  four  possible  wave  configurations 
as  shown  in  Table  1.  R,  C  and  S 


wave  1 

wave  2 

wave  3 

R 

C 

S 

S 

C 

R 

R 

C 

R 

S 

C 

S 

Table  1:  Possible  Wave  Configurations 


indicate  whether  the  wave  is  a  rarefaction,  contact  surface  or  shock.  Two- 
rarefaction  and  two-shock  configurations  can  occur  when  waves  interact  with 
each  other  or  reflect  at  boundaries. 

C.  Important  Wave  Relations 

Recall  that  W^and  W4  are  given  by  Equation  (  7 ).  To  solve  the 
Riemann  problem  we  must  find  the  unknown  states,  W2  and  ,  that  lie 
within  i?* .  Algebraic  relations  are  derived  from  physical  principles  in  any 
fluid  dynamics  text  and  can  be  used  to  connect  known  states  to  the  unknown 
states,  W2  and  ,  across  waves  1  and  3  respectively  (Anderson,  1982; 

Courant  and  Friedrichs,  1948).  E.F.  Toro  has  developed  relations  that 
connect  these  states  via  the  following  particle  velocity  equations  (1989): 

“2  =Wi-/i(;?..^i)  (26) 

and 
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«3  =«4+/4(A»^4)- 


(27) 


Recall  that  the  outer  waves  can  be  either  rarefactions  or  shocks.  This  means 
the  definition  of  the  pressure  functions,  /,  and  /, ,  is  dependent  upon  the  wave 

types.  If  the  outer  wave  is  a  shock,  the  Rankine-Hugoniot  relations  are  used 
to  derive  these  functions.  If  the  wave  is  a  rarefaction,  the  compatibility 


equations,  called  Riemann  invariants,  along  the  C  and  characteristics 
and  isentropic  relations  are  used.  Application  of  these  relations  results  in  the 
following  generalized  pressure  functions  (Toro,  1997): 


f4{p*,W^)  = 


(P*  -P4) 


p^(p*(r+i)+p^(r-i)) 


if  p,  >  p^  shock 


=  < 


2  c. 


y-1 


'Pl" 


Y-l/ 


/2y 


(28) 


if  p*  <  p^  rarefaction 


where  ^  represents  the  known  constant  state,  1  or  4. 

Within  f?* ,  M,  is  constant,  that  is 

«3-«2=0.  (29) 

Substitution  of  Equations  ( 26  )  and  (  27  )  into  Equation  (  29 )  yields  Toro’s 
nonlinear  algebraic  pressure  function  (Toro,  1997): 

=  f,{p,,W,)  +  f^{p,,W,)  +  u,-u^  =0.  (30) 

Given  constant  states  1  and  4,  p,  is  found  numerically  using  a  root  solver. 
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Once  p,  is  found,  the  other  state  variables  follow  directly.  The  particle 
velocity  within  the  star  region  can  be  calculated  using  either  Equation  ( 26  ) 
or  (  27  ).  However,  because  p,  is  a  numerical  approximation,  albeit  a  rather 

accurate  one,  obtained  using  a  root  solver,  u,  is  calculated  as  (Toro,  1997): 

=^(“i+“4)+^(/4-/i)-  (31) 

Unfortunately,  density  within  andi?3  is  dependent  upon  the  wave 

types  of  the  outer  waves  and  therefore  no  simple  relation  exists  to  calculate 
these  variables.  We  must  again  use  a  generahzed  density  function  derived 
from  the  same  relations  used  above  to  derive  Equation  (  28 ).  To  calculate  P2 
and  P3,  the  following  generalized  function  is  used  (Toro,  1997): 


P4 


^p*(r+i)+p^(/-i)^ 

,p*(r-i)+p^  (/+!); 


if  p,  > 


shock 


P^ 


if  Pt  <  p^  rarefaction 


(32) 


where  ^  again  denotes  variables  from  the  known  constant  states. 

D.  Toro’s  Exact  Riemann  Solver  (ERS) 

With  these  relations  now  in  hand,  we  seek  the  exact  solution  to  the 
one-dimensional  Euler  equations  given  the  initial  conditions  in  Equation 
( 18  ).  Unfortunately,  no  closed  form  solution  exists.  We  can  however, 
approximate  the  solution  with  as  high  a  degree  of  accuracy  as  our 
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computational  resources  will  allow.  To  find  the  exact  solution,  we  have 
chosen  Toro’s  exact  Riemann  solver  for  its  simplicity  and  computational 
efficiency  (Toro,  1989). 


The  exact  Riemann  solver  consists  of  three  principle  steps,  as  shown  in 
Figure  5,  beginning  with  the  most  computationally  expensive  portion  of  the 


Figure  5:  Exact  Riemann  Solver  Solution  Procedure 


entire  algorithm,  the  star  step.  In  this  step,  the  variables  within  R.  are 
obtained  given  known  states  Wi  and  using  the  wave  relations  presented 
in  Section  C.  Once  W2  and  1^3  are  known,  the  solution  is  sampled  at  a  user- 

defined  number  of  points.  This  is  called  the  sample  step.  If,  during  the  star 
step,  one  of  the  outer  waves  is  identified  as  a  rarefaction,  the  sample  points 
within  the  expansion  fan  are  determined  using  another  set  of  equations. 

The  Star  Step 

Toro’s  exact  Riemann  solver  begins  with  the  star  step.  Notice  from 
Equations  (  31 )  and  (  32  )  that  u* ,  P2  and  are  all  functions  of  the 
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remaining  star  variable,  p* .  This  means  p,  must  be  found  before  the 
solution  can  be  obtained  fully.  The  pressure  within  the  star  region,  p* ,  is 
approximated  by  finding  the  root  of  Equation  (  30  )  iteratively  using  the 
Newton-Raphson  method.  This  method  takes  the  form  (Burden  and  Faires, 
1993): 

dp* 

where  k  denotes  the  iteration.  It  is  important  to  reiterate  here  the  fact 
that  both  pressure  functions,  /j  and  4 ,  are  defined  differently  depending  upon 
the  nature  of  the  outer-waves.  Therefore,  if  we  wish  to  converge  rapidly  to 
the  solution,  waves  1  and  3  must  be  first  be  characterized  and  then  the 
appropriate  pressure  function  applied. 

Given  f  ^  0  and  a  sufficiently  close  initial  approximation,  , 
Equation  (  33 )  will  converge  quadratically  to  the  root.  The  find  an  adequate 
initial  guess,  Toro’s  adaptive  initial  guess  generator  is  used  (see  Appendix  A). 
The  resulting  approximation  is  considered  adequate  when  the  relative  error 
is  less  than  the  user-defined  error,  e ,  (Toro,  1989): 

<£■  ( 34 ) 
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After  p,  is  found,  the  remaining  star  variables  are  calculated  per 


Section  C  and  the  solution  within  constant  regions  1  through  4  is  fully 
obtained.  Now  the  solution  is  sampled  across  the  computational  domain. 

The  Sample  Step 

In  the  sample  step,  we  sample  the  solution  at  a  user-defined  number  of 
points.  The  solution  at  {x,t)  is  obtained  by  exploiting  the  self-similar  nature 

and  wave  structure  of  the  Riemann  problem.  Recall  that  for  the  solution  to 
be  self-similar,  the  Riemann  problem  domain  must  be  symmetric  about  the 
initial  discontinuity  located  at  jc  =  0  as  shown  in  Figure  6.  S^ ,  Sg  and 

Sg  represent  the  speeds  of  waves  1,  2  and  3  respectively.  In  the  laboratory 


Figure  6:  Riemann  Problem  Local  Reference  Frame 
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reference  frame  the  discontinuity  is  located  at  jc  = 


where  L  represents 


the  physical  length  of  the  problem.  The  physical  domain  is  typically 
[OjLjxfOjT]  where  T  is  the  desired  solution  time.  To  transform  from  the 
physical  to  computational  domain,  the  following  coordinate  transformation 
must  be  performed  to  meet  the  conditions  required  for  similarity: 


(35) 


where  x  andx  represent  position  in  the  local  and  laboratory  reference  frames 
respectively. 

Now  we  take  advantage  of  the  similarity  properties  described  in 
Section  B  to  determine  W(jx,T) .  Here  we  define 


(36) 


A 

where  S  is  the  speed  required  to  move  to  x  at  the  time  T  starting  at 
(0m,0s)  in  the  Riemann  problem  reference  frame. 


Given  S ,  the  solution  is  easily  obtained  by  comparing  it  with  the  wave 


speeds.  If  no  rarefactions  are  present,  the  solution  is: 


mx,T)=\ 


w, 

W, 


if  S<S, 

if  Si  <  S  <  S2 

if  S2  <  S  <  S3 
if  S>S, 


where  the  wave  speeds  are  (Toro,  1991): 


(37) 
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(38) 


Ui  -  Cl 


Y  +  lp,  ^  y-l 

2r  Py  2r 


-il/ 


if  P.  >  Pi 


shock 


> 


-  Cj  if  p,  <  Pi  rarefaction 


and 


'S'2  =  u, , 


(39) 


^3  = 


M4  +C4 


r+1  p.  ,  r-i 

,  2r  P4  2;^ 


-,1/ 


P.  >  P4  s/iocfe 


^4  +  C4  i/  p,  <  P4  rarefaction 


(40) 


If  a  rarefaction  is  present,  the  solution  within  the  expansion  fan  is  calculated 
using  a  different  set  of  equations.  This  occurs  in  the  rarefaction  step. 

The  Rarefaction  Step 

As  stated  earlier,  if  a  rarefaction  is  present  special  care  must  be  taken 
because  the  solution  varies  continuously  between  the  head  and  tail  of  the 
expansion  fan.  Figure  7  depicts  a  typical  solution  of  the  Riemann  problem 
that  consists  of  a  left  traveling  rarefaction,  right  traveling  contact  surface 
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Figure  7:  Typical  Wave  Structure 
and  a  right  traveling  shock.  The  solution  is  desired  at  point  (x,  T) .  We 

determine  it  lies  within  the  expansion  fan  using  the  following  relation: 

ghead  ^  g  ^  gtail 

where 

ahead 

—  ll-^  —  Cj 

and 

1  —  U*  —  C2  . 

To  determine  the  solution  within  a  left  rarefaction,  the  following  set  of 
equations  is  used  (Toro,  1997): 


(41) 


(42) 


(43) 
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Wf°^{S)  = 


p  =  Pi 


^  ^  ^  (Ui -S) 


-,2/ 


u  = 


7+1  (r+i)ci 
r-i 


/^_i 


p  =  Pi 


y  +  1 
2 


Ci  + 


Uj  +  S 


V 


+  :  ^  («i-S) 


r+1  (r+i)ci 


2r/ 


7-1 


(44) 


Similarly,  we  use  the  following  relations  to  determine  if  (jc,7’)lies  within  a 


right  rarefaction  (Toro,  1997). 


Sf <S<S^ 


head 


(45) 


Shead  „ 

3  -  U4  +  C4 


3  -  U*  +  C3 . 


(46) 


(47) 


If  a  right  rarefaction  is  present  then  the  following  set  apphes  (Toro,  1997): 


W^^{S)  = 


P  =  Pa 


u  ■■ 


r-1 


-(U4  -S) 


iX- 


y  +  1  (y+  l)c4 

2  f  r~l  A 

~  C4  H - 1/'4  +0 

V  2 


P=P4 


y  +  1 
2 


r-1 


lr+1  (r+i)c4 


(ua-S) 


2r/ 


7-1 


(48) 


The  entire  solution  of  the  Riemann  problem  is  now  defined. 

E,  The  Sod  Shock  Tuhe  Problem 

To  illustrate  a  t5^ical  solution  to  the  Riemann  problem,  a  mild  shock 
tube  problem  is  solved  using  the  exact  Riemann  solver  with  10001  points 
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(Figures  8-11).  The  solution  consists  of  a  left  rarefaction,  a  contact  and  a 
shock  on  the  right.  The  wave  structure  is  similar  to  that  depicted  in  Figure 
7.  The  shock  tube  is  one  meter  long  with  a  diaphragm  positioned  at 
X  =  0.5m .  The  following  initial  conditions  exist: 


VF(a:,0)  = 


Wi  =il.000kg/m^  ,0.0m/s,1.000PaY  x<0.5m 
-  /  ;  o  vr  ■  ( 49 ) 

1^4  =(0.125fe^/m'*,0.0m/s,0.1OOPaj  x>0.5m 


This  problem  was  first  presented  and  used  by  Sod  (Sod,  1978). 

At  t  =  0s,  the  diaphragm  is  instantaneously  removed  resulting  in 


particle  flow  from  left  to  right  down  the  pressure  gradient.  As  this  flow 


continues  it  accelerates  such  that  particle  velocities  become  supersonic 
relative  to  the  ambient  conditions  or  quiet  conditions  downwind  (Figure  9). 
As  a  result,  compression  waves  coalesce  forming  a  shock  that  travels  to  the 


right  into  the  low-pressure  region  (Figure  10).  Directly  behind  the  shock,  a 


contact  surface  follows  traveling  at  a  lower  speed  (Figure  8).  The  contact 
surface  occurs  because  the  pressure  and  particle  velocity  behind  the  shock  is 


constant  within  i?* .  This  results  in  a  difference  in  internal  energy  between 
the  colder  air  to  the  left  of  the  contact  discontinuity  and  the  shock-heated  air 
to  the  right  as  observed  in  Figure  11.  In  the  opposite  direction,  a  rarefaction 
wave  travels  to  the  left  into  the  high-pressure  region.  Notice  in  Figure  9  that 
the  particle  velocity  is  positive  and  increases  across  the  rarefaction  wave. 

The  expansion  fan  occurs  due  to  the  low  density  region  created  locally  as  the 
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p(kg/m^3) 


heated  gas  moves  rapidly  to  the  right  down  the  pressure  gradient.  This 
results  in  a  region  where  pressure,  density  and  temperature  decrease. 


1.20 


x(in) 

Figure  8:  Density  Profile  at  t=0.25s 
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Figure  11:  Specific  Internal  Energy  Profile  at  t=0.25s 
To  test  the  implementation  of  the  exact  Riemann  solver,  we  solved  five 
shock  tube  problems  in  addition  to  Sod’s: 

-  Toro’s  shock  tube  problem  (1991), 

-  two-rarefaction  shock  tube  problem  (Einfeldt  and  others,  1991), 

-  left  half  of  Woodward  and  Colella’s  blast  problem  (1984), 

-  right  half  or  Woodward  and  Colella’s  blast  problem  (1984)  and 

-  two-shock  shock  tube  problem  (Toro,  1997). 

In  each  case,  we  found  the  results  to  be  in  agreement  with  those  presented  in 
the  literature. 
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Chapter  4:  One-dimensional  Shock  Code  Implementation 

While  the  exact  Riemann  solver  provides  an  exact  solution 
unfortunately,  it  solves  only  the  simplest  of  problems,  namely  symmetric  one¬ 
dimensional  shock  tube  t5T)e  problems.  To  solve  more  interesting  problems 
we  must  turn  to  finite  volume  techniques.  Here  we  describe,  implement  and 
test  several  numerical  techniques  to  build  a  second-order  accurate,  one¬ 
dimensional  shock  code. 

A.  Initial  Boundary  Value  Problem 

In  this  chapter  we  will  consider  and  subsequently  solve  the  following 
general  initial  boundary  value  problem  (IBVP): 


IBVP  = 


Ut+F(U)^=6 
[U(x,0)  =  U^(x) 


(50) 


Unlike  the  RP  described  in  the  last  chapter,  no  symmetry  is  assumed.  In  the 
most  general  sense,  t/®  (jc)  is  a  piecewise  function  that  defines  multiple 
discontinuities  at  t  =  Os  over  the  spatial  domain,  [0,  L] .  The  boundary 
conditions  along  x  =  Om,  and  x  =  L  are  explicitly  defined  to  be  either 
transmissive  or  reflective.  These  boundary  conditions  arise  from  the 
numerical  requirement  to  model  flows  along  physical  and  computational 
boundaries. 
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B.  Discretization  of  the  Domain 

To  solve  Equation  ( 50  ),  the  domain,  [0,L]  x  [OjIT] ,  is  discretized  in  the 
x-t  plane  forming  a  rectilinear  mesh  as  shown  in  Figure  12.  Mesh  indices  i 
and  n  represent  the  cell  at  the  n***  time  level.  The  spatial  domain  is 


t 


Figure  12:  Discretized  Domain 

subdivided  into  I  computational  cells  of  uniform  width  Ax  such  that: 


(51) 


The  center  of  the  computational  cell  is  located  at: 


Xi  =L 


i-0.5 


for  i  =1, 2, .../. 


(52) 


The  extent  of  the  cell  is  defined  by  faces  i-1/2  and  i  +  1/2  located  at 
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and 


^1-1/2  =  ^ 


^i-1^ 


V  i  y 


(53) 


where 


^i+l/2  ~  ^ 


(64) 


^  ^i+1/2  ^t-1/2  •  (  55 ) 

Notice  in  Figure  12  that  each  time  step,  At ,  is  not  constant  but 
determined  adaptively  (Toro,  1997): 

,71  C  Ax 


forn  =  0,l,  ...N. 


0  71 

*^max 


(56) 


where  C  is  called  the  Courant  coefficient  and  is  the  estimated 


maximum  wave  speed  across  the  spatial  domain  at  time^"^ . 

We  now  have  a  computational  mesh  that  is  fully  defined  with  mesh 
points  located  at  the  center  of  each  computational  cell  at  every  time  level.  As 
we  numerically  solve  Equation  ( 50 ),  cell-averaged  values  of  the  solution, 
U(x,t) ,  are  assigned  to  mesh  points  and  are  denoted: 

U^^U{xi,n.  (57) 

C.  Conservative  Discretization  of  the  Euler  Equations 

Recall  from  Chapter  3  that  the  solution  to  the  Euler  equations  is  not 
smooth  ever5rwhere.  We  have  shown  that  discontinuities  occur  at  the  shock 
fi’ont  and  contact  surface.  In  these  regions,  the  solution  is  not  differentiable 
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and  the  differential  form  of  the  Euler  equations  fails.  The  PDEs  are 
integrable  over  these  same  discontinuous  regions  so  we  recast  the  PDEs  into 
the  following  form: 


(58) 

/ace£_i/2  /acej+i/2 


i 

i  1 

1 

1 

1 

1 

1 

\ 

1 

1 

1 

1 

1 

!  ^t+l/2  “ 

1 

1 

-►  j 

1 

r  ! 

1 

1 

I 

1 

1 

1 

1 

1 

Xi- 

-1/2  Xi 

^i-^l/2 

Figure  13:  Arbitrary  Mesh  Rectangle 

If  we  apply  Green’s  theorem  in  the  x-t  plane  and  integrate  Equation  ( 58 ) 
over  an  arbitrary  rectangle,  shown  in  Figure  13,  the  conservative  discretized 
form  is  obtained: 

(59) 

and  Uf  are  space-integrated  averages  of  the  solution  over  cell  i  where 

(60) 

Ax  ■’^.-1/2  ^  ' 

and 

Uf  (61) 

Ax  •’*1-1/2 
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Similarly,  and  >  represent  the  time-integrated  averages  of  the 

physical  flux  of  U  through  cell  faces  i-1/2  and  i  +  1/2  overtime  step  Af . 
These  flux  functions  are  obtained  in  the  following  manner: 

Fi-y2  =  —  F[U(Xi_y2,t))dt  ( 62  ) 

and 

.  ( 63  ) 

Equation  (  59 )  is  a  direct  mathematical  statement  of  the  conservation 
laws.  It  describes  the  variation  in  cell-averaged  variables  over  At  which 
result  from  the  balance  of  the  time-averaged  fluxes  at  the  cell  faces.  The 

solution  procedure  is  obvious,  t/f  is  calculated  from  t/f  and  the  net  flux  , 
(^1-1/2  “-^1+1/2 )» through  cell  i.  Unfortunately,  .Fj_j/2and  .Fi+1/2  are  not 
known  a  priori.  Therefore,  we  need  a  numerical  scheme  that  approximates 
the  flux  through  cell  i  given  Uj^_i  and  Uf . 

D.  Godunov’s  Method 

To  solve  the  IBVP,  we  use  Godunov’s  conservative  upwind  scheme 
(LeVeque,  1992;  Hirsch,  1990).  This  method  assumes  the  solution,  U(x,t), 
has  a  piecewise  constant  distribution  of  cell  averages  over  the  spatial  domain 
at  each  time  level.  At  f  =  Os,  we  map  U^{x)  onto  our  computational  domain 
using  cell-averaged  values  as  defined  by  Equation  ( 61 ).  For  example,  let’s 
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assume  a  discontinuity  exists  at  x  =  as  depicted  in  Figure  14.  Then  the 
cell-averaged  value  assigned  in  cell  3  is  simply  =  (t/J  +  C7j)/2. 


Un. 


Figure  14:  Component  of  U(x,t^)  Mapped  onto  Discretized 

Spatial  Domain 

At  this  point,  we  have  defined  a  set  of  I  constant  states  that  may  be 
interpreted  as  7  - 1  pairs,  |.  Each  of  these  pairs  is  separated  by  a 

discontinuity  that  occurs  at  the  shared  cell  face.  Recall  Fi_]/2  and  Fi+1/2  are 

not  known  but  are  required  if  the  solution  at  is  to  be  found. 

Fortunately,  our  interpretation  of  the  solution  distribution  suggests  a 
simple  procedure.  We  can  calculate  the  unknown  fluxes  by  defining  and 
solving  a  series  of  local  Riemann  problems  exactly  at  each  cell  face: 
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(64) 


IVP  = 


U^+F{U)^=0 

tJixfi) 


where 


t7(x,0)  = 


jc  <  0 
X  >  0 


(65) 


Recall  that,  in  practice,  we  solve  the  Riemann  problem  in  primitive  variable 
form  as  stated  in  Chapter  3.  From  the  exact  solution  of  Equation  ( 64 ),  we 
obtain  the  constant  solution  along  face  i  + 1/2 .  Evaluation  of  Equations  ( 62  ) 
and  ( 63 )  is  now  trivial.  In  general,  the  Godunov  flux  is: 

(66) 

where  lF(Xj^.j/2)  is  dependent  upon  position  only  and  represents  the  primitive 


variable  state  along  face  i  +  1/2 .  Once  Fi_y2  and  Fj+j/2  are  known,  we 
obtain  using  Equation  (  59 ). 

E.  Implementation  of  Godunov’s  Method 

To  calculate  the  solution  at  from  t”' ,  we  implement  Godunov’s 
method  using  six  steps  as  shown  in  Figure  15.  In  step  1,  boundary  conditions 
are  applied  to  fictitious  cells,  0  and  7  + 1 ,  called  phantom  cells.  As  we  shall 
see,  these  cells  are  required  to  obtain  the  solution  in  border  cells  1  and  I. 
Once  boundary  conditions  have  been  applied.  At  is  determined  subject  to  the 
Courant-Friedrichs-Lewy  stability  condition. 
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Figure  15:  Godunov  Solution  Procedure 


In  step  3,  we  solve  all  7  +  1  RPs  to  obtain  along  each 

physical  cell  face.  By  physical  we  mean  those  faces  or  cells  that  lie  within  the 
actual  spatial  domain,  i.e.  no  phantoms.  Most  of  our  computational  effort 
and  time  is  spent  in  this  step  solving  Riemann  problems  exactly.  Step  4 
handles  rarefactions  where  they  occur.  Once  the  solution  is  obtained  along 
each  cell  face,  we  calculate  the  Godunov  flux  in  step  5  using  Equation  ( 66 ). 
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Finally,  we  calculate  t/f  in  each  cell  using  Equation  (  59  ).  Next,  we  will 

look  at  steps  1,  2,  4  and  6  in  greater  detail.  Steps  3  and  5  are 
straightforward. 

Step  1:  Phantom  Cells  and  Boundary  Conditions 

Until  now,  we  have  ignored  the  numerical  difficulties  we  face  along  the 
boundaries  of  our  spatial  domain.  Namely,  border  cells,  1  and  I ,  cannot  be 
updated  using  Equation  ( 59  )  without  Fq  and  .  To  define  these  fluxes,  we 
apply  boundary  conditions  along  x  =  Om  and  x  =  L. 

Traditionally,  these  boundary  conditions  manifest  themselves  in  the 
form  of  boundary  functions  that  explicitly  define  the  unknown  quantities  of 
interest.  Here,  we  take  a  different  approach  and  define  cell-averaged  states 


and 


within  phantom  cells  as  depicted  in  Figure  16  (Toro,  1997). 


(68) 
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Figure  16:  Phantom  Cells 

In  modeling  shocks,  we  consider  only  two  types  of  boundaries:  physical 
boundaries  that  are  reflective  and  computational  boundaries  that  are 
transmissive.  To  model  the  reflection  of  shocks  and  other  waves  at  a 
perfectly  reflective  physical  barrier,  we  set  the  primitive  variable  states 
within  each  phantom  cell  in  the  following  manner: 

1^^  =  -Ui  for  a  reflective  boundary  along  x  =  Om  ( 69 ) 

kPi  > 

and 

Pi 

-Uj 

Pi 
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The  solution  at  reflective  boundaries  will  consist  of  either  two-shocks  or  two- 
rarefactions. 

Transmissive  boundaries  are  applied  to  restrict  our  physical  problem 
to  a  finite  computational  domain.  Ideally,  the  flows  should  propagate 
through  hypothetical  boundaries  without  any  effect.  For  transmissive 
boundaries  we  simply  set  the  unknown  phantom  state  equal  to  that  in  the 
corresponding  border  cell,  i.e. 

Wq  =  for  a  transmissive  boundary  at  x  =  Om  (  71 ) 

or 

^j+i  =  for  a  transmissive  boundary  at  x  =  L .  (  72  ) 

This  technique  works  reasonably  well  in  one-dimension.  In  multiple 
dimensions  however,  noticeable  computational,  or  false,  reflections  occur. 

The  inability  to  model  transmissive  BCs  in  multiple  dimensions  has  led  to  a 
significant  and  ongoing  research  effort. 

With  the  state  of  the  phantom  cells  determined,  local  Riemann 
problems  may  now  be  solved  along  the  boundaries  to  calculate  the  unknown 
fluxes,  Fq  and  . 

Step  2:  Determination  of  a  Stable  Time  Step 

Before  we  can  move  forward  in  time  to  solve  Riemann  problems  and 
update  the  solution  an  appropriate  time  step  must  be  determined  that  meets 
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the  Courant-Friedrichs-Lewy  stability  requirement  (Courant  and  Friedrichs, 
1948): 

(73) 

Ax 

where  <S^ax  the  maximum  wave  speed  that  occurs  across  the 

computational  spatial  domain.  The  importance  of  this  condition  in  properly 
selecting  At  becomes  evident  if  we  consider  adjacent  RPs.  In  Figure  17,  the 
solution  to  the  left  Riemann  problem  has  a  strong  shock  wave  traveling 
rapidly  to  the  right.  If  At  is  too  large  as  depicted  in  Figure  17,  the  shock 
wave  travels  far  enough  to  the  right  that  it  perturbs  the  solution  we  seek 
along  face  i  +  3/2 .  This  violates  the  assertion  made  earlier  that 


Figure  17:  Wave  Propagation  Through  Adjacent  RPs 
W(xi+y2)  is  constant.  This  will  result  in  an  incorrect  flux  calculation 

creating  numerical  errors  that  propagate  outward  affecting  other  cells.  As  a 
result,  the  code  rapidly  becomes  unstable  and  crashes.  To  prevent  this,  the 
time  step  is  reduced  to  At  as  shown. 


44 


Ideally,  we  would  like  to  maximize  the  computational  efficiency  of  the 
code  by  taking  the  largest  time  step  possible: 


At  = 


Ax 


qrZl 

^max 


(74) 


where  satisfies  Equation  (  73  )  exactly.  Unfortunately,  the  maximum 
wave  speed  is  not  known  a  priori  but  must  be  estimated  (Toro,  1997): 

‘^max  ~  ^^Sx(|ltj|  +  Cj)  .  (  75  ) 

This  estimate  is  chosen  because  it  extends  easily  into  multiple  dimensions. 
However,  it  requires  some  caution. 

During  the  first  few  time  steps,  our  estimate  is  dominated  by  the  local 
speed  of  sound  and  results  in  a  significant  underestimate  of  S^ax  •  If  we 

simply  substitute  iiito  Equation  (  74 )  we  will  calculate  a  At  that  is  too 

large  and  the  code  will  become  unstable.  To  prevent  such  an  occurrence,  we 
introduce  a  scaling  factor,  C ,  called  the  Courant  coefficient  into  Equation 
(  74 )  that  satisfies  0  <  C  <  1 .  Each  time  step  is  then  calculated  adaptively 
using  Equation  (  56 ).  During  the  first  several  time  steps,  Toro  suggests 
C  =  0.2  (1997).  At  later  time  steps,  the  Courant  coefficient  is  set  to  0.9  to 
increase  At  and  computational  efficiency  while  maintaining  stability. 

Step  4:  Sonic  Rarefactions  -  A  Special  Case 

When  solving  RPs,  we  must  consider  a  special  case  -  sonic  rarefactions. 
These  rarefactions  occur  when  the  head  and  tail  of  the  expansion  fan  travel 
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in  opposite  directions.  When  this  occurs,  the  solution  lies  within  the 
expansion  fan.  In  order  to  calculate  the  solution,  we  must  determine  the 
directions  of  travel  and  apply  the  appropriate  rarefaction  equation.  Equation 
( 44 )  or  (  48 ),  presented  in  Chapter  3. 

Step  6:  Update  Cell  Solution 

At  this  point,  the  solution  to  all  the  Riemann  problems  are  known,  the 
solution  obtained  along  each  cell  face  and  the  corresponding  fluxes 
calculated.  It  is  important  to  remember  that  while  the  primitive  variable 
form  is  used  in  solving  the  Riemann  problems  and  calculating  the  unknown 
fluxes,  only  the  conservative  formulation  will  result  in  an  accurate  solution. 

The  solution,  t/f  ,  is  calculated  using  Equation  ( 59 )  given  and  fluxes 
Fj_]/2  and  .^+1/2  •  Before  the  next  time  step  begins,  is  calculated  from 

0n+l 

F.  Illustration  of  Godunov’s  Method 

To  illustrate  the  important  numerical  properties  of  Godunov’s  method, 
we  revisit  the  classic  shock  tube  problem  presented  in  Chapter  3.  Recall  that 
the  solution  consists  of  a  left  traveling  rarefaction,  right  traveling  contact 
surface  and  a  right  travehng  shock  wave.  The  solution  was  obtained  after 
125  time  steps,  t  =  0.25s ,  using  a  course  computational  mesh  of  200  cells. 
During  the  first  five  time  steps,  the  Courant  coefficient  was  set  to  0.2.  After 
that,  it  was  set  to  0.9.  Boundary  conditions  were  transmissive. 
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We  present  the  density  profile  obtained  via  Godunov’s  method  in 
Figure  18.  For  comparison,  the  exact  solution  obtained  in  Chapter  3  is  shown 
as  well.  The  solution  is  smeared  near  each  discontinuity.  This  undesirable 


Figure  18:  Godunov  vs  ERS  -  Density  Profile  at  t=0.25s 
behavior  is  common  among  first-order  methods  and  is  caused  by  excessive 
numerical  viscosity. 

Clearly,  this  smearing  is  not  equal  in  severity  everywhere.  Notice  that 
the  shock  front  is  spread  over  six  cells  while  the  contact  surface  occurs  over 
23  cells.  This  disparity  in  resolution  is  common  among  Riemann  problem 
type  methods  and  is  a  direct  result  of  the  characteristic  nature  of  each  wave. 
For  a  contact  wave,  the  characteristics  run  parallel  to  one  another  and  it  is 
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the  difference  in  internal  energy  that  results  in  the  density  discontinuity. 
Therefore,  the  contact  discontinuity  is  not  easily  resolved.  Shock  wave 
resolution  is  better  because  the  characteristics  run  into  each  other.  The 
smearing  that  occurs  along  the  rarefaction  fan  occurs  at  the  head  and  tail 
where  the  first  derivative  of  the  solution  is  discontinuous.  The  rarefaction 
wave  becomes  resolved  quickly  as  we  refine  the  mesh. 

Before  we  leave  this  section,  several  positive  features  of  Godunov’s 
method  are  pointed  out.  The  position  and  speed  of  each  wave  are  predicted 
accurately.  This  is  very  important  because  we  want  to  model  shocks 
correctly.  The  monotonic  behavior  of  the  solution  near  discontinuities  is 
another  important  quality.  This  feature  becomes  important  later  when  we 
apply  higher-order  techniques.  Lastly,  while  numerical  viscosity  is  evident, 
most  first-order  methods  result  in  significantly  more  smearing  than  we 
observe  here. 

G.  Improvements  in  Computational  Efficiency  and  Accuracy 

Adaptive  Riemann  Solver  (ARS) 

In  Section  F,  we  found  Godunov’s  method  to  be  first-order  accurate. 
Given  the  large  computational  effort  spent  solving  each  Riemann  problem 
exactly,  this  method  seems  dubious  at  best.  The  obvious  question  comes  to 
mind.  Why  put  forth  all  that  effort  to  achieve  only  first-order  accuracy? 
Fortunately,  over  the  last  two  decades,  several  algorithms  have  been 
developed  that  approximate  the  solution  to  the  Riemann  problem  with 
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accuracy  and  efficiency.  These  approximate  Riemann  solvers  have  given  rise 
to  an  entire  class  of  techniques  called  Godunov-type  methods.  Approximate 
Riemann  solvers  are  used  in  one  of  two  ways: 

-  they  are  used  to  approximate  the  numerical  flux  directly 

-  they  are  used  to  approximate  the  state,  W{Xi+y2)- 

Here,  we  use  an  adaptive  scheme  that  switches  between  three  approximate 
Riemann  solvers  to  find  a  state  approximation.  We  then  calculate  the 
Godunov  flux  and  proceed  as  described  in  Section  E. 

The  adaptive  Riemann  solver  scheme  was  first  introduced  by  Toro  in 
1997.  We  have  implemented  this  particular  algorithm  because  it 

is  easily  implemented  within  our  code  —  two  of  the  approximate 
Riemann  solvers  are  derived  directly  from  the  exact  Riemann 
solver, 

-  is  very  efficient  computationally  —  none  of  the  approximate 
Riemann  solvers  require  iteration, 

-  captures  shocks  well  -  does  not  require  shock  fitting  or  entropy  fix 
common  among  other  techniques  and  it 

-  predicts  the  appropriate  solver  for  each  RP  adaptively  using  a 
technique  we  have  already  employed  and  found  successful  (see 
appendix  A). 
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As  stated  above,  ARS  selects  one  of  three  approximate  Riemann  solvers  based 
upon  the  local  flow  conditions  present  and  the  behavior  of  Toro’s  pressure 
function.  In  fact,  ARS  uses  the  same  switching  criteria  introduced  in 
appendix  A  to  predict  an  initial  guess  for  the  Newton-Raphson  root  solver. 

We  summarize  that  logic  here: 

1.  Calculate  and  Pp^^.^. 

2.  If  Pmiu  —  Ppvrs  —  Pmax  and  Qg  —  Pmax/ Pmin  then  lies  within  1 2 
and  mild  conditions  exist.  Approximate  using  the  primitive 
variable  solver. 

3.  If  however,  Pp^^s  <Pmin  P*  lies  within  and  two  rarefactions 
are  present.  Approocimate  i?,  using  the  two-rarefaction  solver. 

4.  Else,  p,  lies  within  J3  and  either  two  shocks  or  a  single  strong  shock 
are  present.  Approodmate  R*  using  the  two-shock  solver. 

The  three  approximate  Riemann  solvers  of  the  adaptive  Riemann  solver  are: 

-  the  primitive  variable  Riemann  solver, 

-  the  two-rarefaction  Riemann  solver  and 

-  the  two-shock  Riemann  solver. 

We  present  the  solution  method  below  for  each. 

Primitive  Variable  Riemann  Solver  (PVRS) 

In  typical  shock  problems,  most  of  the  flow  field  varies  smoothly.  If  the 
solution  lies  within  I2  and  no  strong  shocks  are  present,  we  may 
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approximate  the  solution  by  solving  a  linear  hyperbolic  system  of  PDEs 
exactly.  Toro’s  primitive  variable  solver  is  derived  in  this  fashion  and  results 
in  a  set  of  simple  algebraic  expressions  (1991).  Using  PVRS,  we  approximate 
the  solution  in  i?*  using  the  following  relations  (Toro,  1991): 


P*  +P4)  +  ^(“i  -M4)(A  +  A)(Ci  +C4)  , 


2  (/?]  +  +C4) 


P2  -  Pi  + 


(«i  -»>)(Pi  +P4) 
(c,  +C4) 


(76) 


(77) 


(78) 


and 


Ps  =P4  + 


(«. -M4)(Pi 

(c,  +C4) 


(79) 


Two-Rarefaction  Riemann  Solver  (TRRS) 

If  two  rarefactions  are  present  within  the  solution,  we  can  find  p, 
exactly  using  the  appropriate  rarefaction  relations  in  Equation  ( 30 ).  Using 
the  two-rarefaction  solution,  we  approximate  the  solution  within  JR,  as  (Toro, 


1997): 


Cl  +C4  --(;'-1)(M4  -Uj) 

_  ^  _ 


(y-i)/2r  „  (x-i)/2r 

V  Pi  P4 


(80) 
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(81) 


M.  =^(M,  +M4)+|(/4-/i) 


and 


d^±iiP  *  i  P^iP^)  - 


P4 


^p*(r  +  i)  +  p^(/-i)^ 
p*(r-i)+p^(/+i) 


P^ 


/  \ 
p* 


\Po 


if  p,  >  shock 


if  p*  <  p^  rarefaction 


(82) 


Two-Shock  Riemann  Solver  (TSRS) 

If  two  shocks  or  a  single  strong  shock  are  present,  we  use  the  two- 
shock  approximation.  Direct  substitution  of  the  appropriate  shock  relations 
into  Equation  (  30 )  results  in  the  two-shock  approximation.  Unfortunately, 
this  approximation  is  dependent  upon  p,  itself  and  therefore,  does  not  have  a 
closed-form  solution.  An  acceptable  solution  however,  is  to  estimate  p*  first 
using  the  primitive  variable  result  described  above.  This  requires  no 
additional  effort  since  we  have  already  performed  this  estimate  in 
determining  which  solver  to  use.  We  simply  substitute  this  estimate  into  our 
two-shock  approximation.  We  now  have  an  adequate  approximation  without 
iteration. 

We  approximate  the  variables  using  the  two-shock  solver  as  (Toro, 

1997): 
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Pi 


Pi(pl®^(r+i)+Pi(r-i). 


Pa 


N  V 


P4(p*‘^^(r+i)+P4(r-i), 


+  Ui  -U4 


(83) 


Pi(p*®\r+i)+Pi(x-i) 


^r2  / 
+1 


-1-1 


p4(p»®^(r+i)+P4(?'-i), 


where 

=  max(l  X 10"®  Pa,  )  ( 84 ) 

To  obtain  the  most  accurate  solution,  we  calculate  the  remaining  star 
variables  as  described  above  in  TRRS. 

Adaptive  Riemann  Solver  Results 

We  return  to  Sod’s  problem  and  solve  it  again,  using  the  initial 
conditions  and  choice  of  C .  In  Figure  19,  we  compare  the  density  profile 
obtained  using  Godunov’s  method  with  both  the  exact  and  adaptive  solvers. 
Graphically,  we  can  not  differentiate  between  the  two.  Each  results  in  the 
same  first-order  behavior  described  in  Section  F. 
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r  X  Godunov-ARS  i 


X  (m) 

Figure  19:  Comparison  of  Godunov's  Method  Using  ERS  and  ARS 

Computationally  however,  the  difference  is  rather  significant. 
Choosing  the  adaptive  routine  results  in  a  33%  increase  in  computational 
efficiency  as  measured  by  the  amount  of  CPU  time  spent  in  solving  the 
problem,  shown  in  Table  2.  Note  that  many  of  the  Riemann  problems  are 
solved  using  the  primitive  variable  method.  This  is  the  simplest  of  the  three 
approximate  solvers  used  in  the  adaptive  scheme.  Additional  test  cases  and 
problems  have  been  solved  using  both  solvers  and  the  results  observed  are 
consistent  in  every  case.  From  these  results,  we  conclude  the  adaptive  solver 
outperforms  the  exact  one  and  will  be  used  as  the  Riemann  solver. 
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Method 

ERS 

PVRS 

TRRS 

TSRS 

- - 

% 

% 

% 

% 

ERS 

1.28 

100 

0 

0 

0 

ARS 

0.859 

0 

90 

5 

5 

Table  2:  Comparison  Between  the  ERS  and  ARS 


Adaptive  Courant  Coefficient 

In  an  effort  to  improve  computational  efficiency  and  accuracy,  we 
briefly  investigate  an  adaptive  routine  that  sets  the  Courant  coefficient,  C, 
adaptively.  The  original  goal  of  this  effort  was  to  maximize  the  time  step 
from  the  very  beginning  while  maintaining  numerical  stabihty  over  the 
entire  temporal  domain.  Traditionally,  C  is  set  to  a  small  positive  number, 

often  C=  0.2,  to  correct  the  underestimate  of  •  After  several  time  steps, 
usually  5,  C  is  set  to  0.9  to  maximize  our  time  steps  while  maintaining 
stability.  Our  concern  in  using  the  simplistic  approach  above  is  twofold: 

why  waste  computational  efficiency  in  the  early  stages  of  code 
execution  -  can  we  increase  C  sooner,  say  after  one  time  step? 
once  we  set  C  to  0.9  or  some  other  value  determined  adaptively, 
can  complex  flows  ever  result  in  a  loss  of  stability? 

Here,  we  investigate  two  possible  adaptive  techniques  to  maximize  the 
time  step  while  maintaining  stabihty.  At  the  beginning  of  each  time  step,  we 
assume  the  value  of  C  scales  the  right  hand  side  exactly  resulting  in  the 
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maximum  allowable  time  step.  Once  all  the  Riemann  problems  are  solved, 
we  calculate  the  maximum  local  Courant  number, 


^max  =  max 


Ax 


(85) 


If  ^max  1,  the  time  step  satisfies  Equation  (  73  ),  we  reset  C  and 
continue  our  calculations.  It  is  the  method  chosen  in  resetting  C  that  defines 
the  different  adaptive  methods.  In  the  naive  approach,  we  assume  flows  are 
sufficiently  established  over  the  current  time  step,  f"to  This  implies 
that,  at  S^ax  ^  reasonable  estimate.  In  this  case,  we  simply  set 

=  Cjnax  where  C^ax  i®  fh®  user-defined  maximum  allowable  Courant 


coefficient. 

A  more  sophisticated  method  assumes  flow  conditions  will  remain 

essentially  the  same  over  the  next  time  step.  This  suggests  s  SlLti  and 

allows  us  to  usG  tliG  currently  obtained  solution  to  scale  C  in  the  following 
manner: 


=  min 


r  - 

^max  ^max  .  n  ^ 
-hVmax?^] 


max 


We  choose  the  minimum  to  handle  the  special  case  where 

max  max 

Given  typical  flow  conditions,  we  expect  C  will  rapidly  approach 
asymptotically. 


(86) 

<1. 
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If  v^ax  >  1  however,  our  assumptions  concerning  C  were 


inappropriate.  In  this  case  we  throw  away  the  solution,  reduce  and  begin 
again.  In  each  adaptive  routine  described  above,  we  reset  C  in  the  following 
manner: 


= 


f 

\  ^max  / 


(87) 


Adaptive  Courant  Results 

Here,  we  solve  Sod’s  problem  via  Godunov’s  method  (with  ARS)  using 
the  typical  method  as  suggested  by  Toro,  the  naive  adaptive  Courant  method 
and  the  sophisticated  adaptive  Courant  method.  We  discretized  the  spatial 
domain  over  800  cells  to  increase  the  number  of  time  steps. 


Method 

CPU  Time 

Time  Steps 

Restarts 

Toro 

5.2075 

490 

0 

naive 

5.218 

486 

1 

sophisticated 

5.238 

487 

0 

Table  3:  Adaptive  Courant  Results  -  Sod  Test 


We  see  in  Table  3  that  Toro’s  and  our  naive  method  3  calculate  stable 
time  steps  throughout  the  temporal  domain.  The  naive  approach  however, 
required  a  restart.  This  occurred  at  the  second  time  step  because  flows  were 

not  yet  established  sufficiently  to  predict  S^ax  resulting  in  an  unstable  time 

step  during  our  first  attempt.  CPU  times  increase  as  we  increase  the 
complexity  of  the  time  step  calculations.  When  compared  to  the  number  of 
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time  steps  saved,  our  results  indicate  that  adaptive  time  step  methods  are 
counterproductive . 

To  investigate  the  behavior  of  C  in  the  presence  of  complicated  flows, 
we  solve  the  classic  Woodward-Colella  blast  problem  that  consists  of  two 
strong  shocks  traveling  in  air  towards  each  other  (Woodward  and  Colella, 
1984).  These  shocks  interact  with  each  other  and  reflect  from  reflective 
boundaries  resulting  in  stagnation,  reflection  and  transmission  of  shock 
waves.  The  problem  begins  at  t  =  Os  with  three  initial  states  at  rest.  These 
states  are  separated  by  two  discontinuities  located  at  :k  =  0.1/n  and  x  =  0.9/n . 
The  physical  domain  is  Im  in  length  and  lies  between  two  reflective  barriers. 
The  domain  is  discretized  over  1000  cells.  The  following  initial  conditions 
are: 


W(x,0)  =  \ 


W  = 

w  = 
w  = 


il.OOkg/m^ 

lil.OOkglm^ 

i^.OOkgf 


0.00  ml  s ,  lOOO.OOPa]^ 
0.00  m/s,  O.OlPaJ^ 
0.00  m/s,  100. OOPaJ^ 


0.0m  <x  <  0.1m 
0.1m  <  X  <  0.9m . 
0.9m  <x<  1.0m 


(88) 


After  0.05s,  the  shocks  have  reflected  from  the  reflective  boundaries  and  we 
terminate  program  execution.  Table  4  shows  our  results  for  the  same  three 
methods  we  considered  above.  Our  results  are  similar  with  those  obtained  in 
Table  3  despite  the  presence  of  complex  flows;  CPU  time  increases  as  each 
case  increases  in  complexity.  Notice  that  the  naive  method  requires  no 
restart.  The  absence  of  a  restart  is  explained  when  we  consider  the  large 
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pressure  gradients  involved  at  early  times.  This  results  in  a  satisfactory 
^max  aftei*  only  one  time  step.  Our  results  suggest  there  is  no  need  for  an 

adaptive  Courant  coefficient  in  the  presence  of  mild  or  strong  flows.  In  fact, 
using  an  adaptive  method  decreases  the  computational  efficiency  of  the  code 
therefore  we  omit  it. 


Method 

CPU  Time 

Time  Steps 

Restarts 

Toro 

40.59 

2155 

0 

naive 

40.81 

2151 

0 

sophisticated 

42.51 

2152 

0 

Tal 

Ijle  4:  Adaptive  Courant  Results  -  Sod  Test 

H.  Toro’s  Weighted  Average  Flux  (WAF)  Method 

Now  that  we  have  implemented  an  efficient  form  of  Godunov’s  method, 
we  present  Toro’s  second-order  extension  called  WAF.  While  WAF  can  be 
applied  equally  to  many  first-order  techniques,  we  apply  it  here  with 
Godunov’s  method  to  achieve  second-order  accuracy.  Its  implementation  is 
numerically  simple  and  straightforward. 

Recall  that  in  Godunov’s  method  we  calculate  the  numerical  flux  using 
the  solution  of  the  Riemann  problem  obtained  along  the  cell  face  only.  Much 
of  the  characteristic  information  obtained  in  solving  the  Riemann  problem 
has  not  been  utilized,  yet.  WAF  however,  approximates  the  total  flux 
through  each  cell  face  using  a  spatial  quadrature  scheme  across  the  entire 
wave  structure  of  the  Riemann  problem.  Using  all  the  characteristic 
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information  allows  us  to  calculate  a  numerical  flux  that  is  more  accurate. 
Empirical  results  show  this  scheme  is  second-order  accurate. 

We  calculate  WAF  in  the  following  manner  (Toro,  1989): 

^•+1/2  =  (  89  ) 

k  '  ^ 

where  Wj^  is  a  numerical  weight  that  represents  the  geometric  extent  of 
region  k  and 

Fk=F(Wk)  (90) 

is  the  flux  through  region  k .  Figure  20  illustrates  WAF  in  the  x-t  plane  for 
the  Riemann  problem  at  face  i  +  1/2 . 


t 


Figure  20:  Weights  of  WAF 

At  we  see  the  Riemann  problem  spatial  domain  is  subdivided 

into  four  subintervals  of  length  Ax .  Notice  that  we  depict  the  outer  waves 

with  only  one  ray  instead  of  two.  Empirical  evidence  has  shown  that  the 
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role  each  partial  flux  plays  in  determining  is  dependent  upon  .  In 

practice,  we  implement  WAF  in  the  following  manner: 

1.  Calculate  weighted  average  state  along  face  i  +  1/2  where 

k 
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2.  Calculate  weighted  average  flxix  from  our  weighted  average  state 


i  e  =  p(  x 

^•‘^•’•^1+1/2  1+1/2''- 


This  form  of  WAF  is  much  easier  to  implement. 

I.  Illustration  of  WAF 

To  illustrate  WAF  we  solve  Sod’s  problem  and  present  the  density 
profile  again.  In  Figure  21,  we  see  improved  resolution  of  discontinuities  over 
that  achieved  using  Godunov’s  method  alone.  It  is  immediately  obvious 
however,  that  oscillation  occurs  near  each  discontinuity.  This  behavior  is 
quite  t5rpical  of  second-order  techniques.  To  explain  the  origin  of  these 
oscillations,  we  return  to  the  definition  of  WAF. 


Figure  21:  WAF  -  Density  Profile  at  t=0.25s 


62 


Close  inspection  of  Equation  ( 89 )  provides  some  insight  into  WAF  and 
the  oscillatory  nature  of  the  solution.  Let’s  consider  the  Riemann  problem  in 
Figure  20.  Notice  that  the  local  t  axis  lies  within  the  second  region  and 

therefore,  F2  =  .  This  illustrates  an  important  fact;  Godunov’s  flux  is 

always  represented  by  one  of  the  fluxes  in  Equation  ( 89  ).  Note  that  it  is  the 
local  Riemann  problem  structure  that  determines  which  flux  corresponds  to 
Godunov’s  flux. 

F  represents  the  upwind  bias  of  WAF  and  is  responsible  for 
stability.  The  other  fluxes  represent  downwind  terms  and  are  responsible  for 
the  increased  accuracy  we  observe  with  WAF.  It  is  these  downwind  terms 
that  cause  the  oscillations  we  observe  near  steep  gradients.  In  the  next 
section,  we  present  an  adaptive  technique  that  controls  the  contribution  of 
these  downwind  terms  near  discontinuities  resulting  in  near  monotonic 
conditions. 

J.  Limited-WAF  (LWAF) 

Here,  we  present  total  variation  diminishing  (TVD)  weight  limiter 
functions  that  restrict  the  downwind  partial  flux  contributions  near 
discontinuities  by  introducing  dissipative  viscosity  (Toro,  1997).  This  results 
in  near  monotonic  conditions  throughout  the  domain. 
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Before  we  continue,  we  briefly  define  what  TVD  means.  A  numerical 
scheme  is  considered  to  be  TVD  in  nature  when  the  total  variation  of  the 
solution  diminishes  over  time  (LeVeque,  1992).  In  mathematical  terms, 

TV{U^^^)<TV{U^)  (96) 

where  the  total  variation  function,  TV ,  of  the  solution  at  <"is  given  by: 

7V(?7")  =  2; 

i 

TVD  numerical  schemes  have  several  important  features: 

-  they  are  TV  stable  meaning  they  are  guaranteed  to  converge  to  the 
appropriate  solution. 

the  numerical  scheme  decreases  toward  a  monotonic  solution. 
Monotonic-like  second-order  numerical  schemes  are  derived  from  these 
principles. 

To  achieve  monotonic-like  conditions  near  discontinuities,  Toro  has 
developed  TVD  weight  limiter  functions,  ^(r,  v) ,  that  are  based  upon  TVD 

flux  limiter  functions,  .  Any  weight  limiter  may  be  derived  given  the 
corresponding  flux  limiter  using  the  following  relation  (Toro,  1997): 

(^(r,  v)  =  1  -  (l  -  |v|)  y/(r)  <7(v)  ( 93 ) 

where  r  is  a  local  flow  parameter  we  will  discuss  later  and  a  =  SGN(v) . 

Here,  we  construct  the  following  four  weight  limiter  functions:  Super  A  (^5^), 
van  Leer  A  (<%^),  van  Albada  A  (^aa)  and  Min  A  Each  is 


TT^ 

^i+1 


■C/f 


(97) 
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constructed  using  Equation  (  98 )  from  the  corresponding  flux  limiter 
functions  (Toro,  1997): 


0 

r  <  0 

2r 

0<r<0.5 

Superbee:  y/g^  =< 

1 

0.5  <  r  <  1 , 

r 

l<r<2 

2 

r>2 

van  Leer:  y/y^  = 

0 

<  2r 

.1  +  r 

r<0 

r>0’ 

0 


van  Alabada:  y/yj^  =  < 


r(l  +  r) 
l  +  r^ 


r  <0 
r>0’ 


(99) 


(100) 


(101) 


and 


Minbee:  y/MB  = 


r<0 
0<r  <1. 
r>l 


(102) 


Before  we  describe  the  weight  hmiters,  it  is  helpful  to  investigate  how  the 
flux  limiters  work. 

Each  flux  limiter  is  calculated  based  upon  the  local  flow  parameter,  r , 
defined  as  (Sweby,  1984): 


r  = 


^Pupwind 


^Plocal 


(103) 
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where  l^Piocai  is  the  density  change  from  right  to  left  across  a  specific  wave  in 
the  Riemann  problem  of  interest  and  t^Pupwind  is  the  change  in  density  across 

the  same  wave  in  the  adjacent  upwind  Riemann  problem.  In  WAF,  we 
calculate  local  flow  parameters  for  each  wave  in  the  Riemann  problem  at  face 
i  +  1/2  such  that 


^1+1/2, fe 


^A  +  l/2-(7,fe 

^A+l/2,fe 


(104) 


where  the  local  density  change  across  each  wave  is  given  by: 

Vi+i/2,fe  =  A+]/2,fe+i  -Pi+\l2,k  •  Figure  22  illustrates  how  ri+xj2,k  is  calculated. 
In  case  1, 

Case  1 :  >  0 


Case  2 :  <  0 

wavei^yz.k  wavei^3i2^k 
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v;^  >  0  meaning  wave’s  are  traveling  to  the  right.  In  this  case,  <t  =  1 


AyC7j_i/2 

identifying  wavei_ii2  k  as  the  upwind  wave  and  rj+1/2  k  = - In  case  2 

’  ^Pi+ll2,k 


however,  v;^  <  0  and  <7  =  -1  so  that  wavei+^i2^k  is  the  upwind  wave  and 


^Pi+3j2  k 

n+i/2,k  =- — ^  - .  In  each  case,  r  is  a  measure  of  how  rapidly  the  density 
’  ^Pi+l/2,k 

gradient  changes  locally.  This  is  a  measure  of  the  smoothness  of  the  flow,  i.e. 
as  r  ->  1 ,  the  flow  becomes  smooth.  In  these  regions  we  include  the 
downwind  flux  contributions  to  achieve  second-order  accuracy.  As 

»  0 ,  the  flow  becomes  discontinuous.  In  these  regions  we  severely 

limit  the  contribution  of  all  second-order  downwind  terms  to  to 


achieve  stabihty  and  monotonic  behavior. 

Weight  limiters  work  in  a  similar  manner.  Instead  of  limiting  the  flux 
contributions  directly  however,  weight  limiters  amplify  courant  numbers,  , 

depending  upon  local  flow  conditions  as  follows: 


<l>k=^k=  OCk  Vk 


(105) 


Kl 


is  the  TVD  wave  amplifier  for  the  wave 


(Toro,  1992).  Amplification  of  these  Courant  numbers  changes  the  value  of 
the  corresponding  weights.  For  example,  as  |l  -r^l  becomes  large,  the 
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corresponding  weight  approaches  zero  effectively  killing  off  the  downwind 
contribution  of  that  flux. 

Boundary  conditions  in  the  limited  WAF  scheme  must  be  treated  with 
care.  In  order  to  account  for  the  upwind  calculations  required  by  the  limited 
WAF  scheme,  phantom  cells  -1  and  1  +  2  must  be  added.  For  example, 
consider  the  calculation  of  the  limited  flux  at  x  =  Om .  If  Vii2,k  >  0  >  then  we 

must  add  cell  - 1  to  solve  the  Riemann  problem  located  at  x  =  -Ax . 

K.  Illustration  of  Limited  WAF 

Now  that  we  have  developed  the  limited  version  of  WAF,  we  solve 
Sod’s  problem  using  each  weight  limiter  function  to  briefly  investigate  the 
numerical  properties  of  each.  The  results  are  shown  in  Figures  23  -  25.  In 
Figure  23  we  see  that  every  limiter  dampens  out  the  second-order  partial  flux 
terms  near  discontinuities.  It  is  however,  difficult  to  pick  out  any  details 
concerning  the  performance  of  each  limiter  near  areas  of  interest:  contact 
discontinuities  and  shocks.  Figures  24  and  25  amphfy  these  regions.  It 
becomes  evident  from  Figure  24  that  Super  A  performs  best  near  contact 
surfaces.  In  Figure  25  the  behavior  of  each  appears  to  be  similar  until  we 
observe  the  behavior  of  each  to  the  left  of  the  shock.  In  this  region  we  see 
that  van  Leer  A  performs  best.  The  other  limiter  functions  oscillate  slightly. 
From  these  results  we  conclude  that  van  Leer  A  is  appropriate  when  we 
model  shocks. 
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Figure  25:  Comparison  of  Limiters  near  Shock 
L.  Testing  and  Validation  of  the  One-Dimensional  Shock  Code 
Testing 

To  test  the  implementation  of  the  one-dimensional  code  we  solved  the 
five  shock  tube  problems  presented  in  Chapter  3  in  addition  to  Sod’s  as 
presented  here.  At  each  step  in  the  development  of  the  code,  we  qualitatively 
compared  our  results  with  those  obtained  using  the  exact  Riemann  solver. 
This  was  done  to  ensure  proper  implementation  of  each  numerical  tool  by 
testing  it  over  a  broad  range  of  problems  before  adding  additional  complexity. 
Having  tested  for  proper  implementation,  we  then  validated  the  code  against 
experimental  data. 
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Validation 


In  1996,  the  Army  Research  Laboratory  used  measurements  obtained 
in  their  51cm  shock  tube  facility  to  validated  their  own  in-house  second-order 
accurate  one-dimensional  shock  code.  Here,  we  use  the  same  data  to  vahdate 
our  own  code  before  moving  into  two  dimensions.  The  shock  tube  was  100m 
in  length,  51cm  in  diameter  and  contained  a  driver  region  0.91m  long 
(Schraml,1996).  A  measurement  station  was  located  31.48m  from  the  initial 
discontinuity  where  overpressure  time  histories  were  taken  experimentally 
for  validation  of  the  Army  Research  Laboratory  code.  The  following 
conditions  were  present  at  the  beginning  of  the  experiment: 

=  (4.486% /m^,  O.Om/s,  379.2 X 10^ Pa)  X  <  0.91m 
W(x,0)  =  \  .  .  (106) 

=  (l.208%/m^,0.0m/s,102.1xl0^Paj  x>  0.91m 

Figure  26  shows  the  overpressure  time  history  taken  at  the 
measurement  station  as  compared  to  the  Army  Research  Laboratory’s 
computational  results  (Schraml,1996).  Notice  the  large  oscillations  that 
occur  near  the  shock  front  in  Figure  26.  The  shock  front  shock  arrived  at  the 
measurement  station  after  66.0ms.  The  measured  peak  overpressure  at  the 
shock  front  was  66.3/ePa. 


72 


60  70  80  90  100 

Time  (ms) 


Figure  26:  ARL  Overpressure  Time  History  at  31.44  m 
To  correctly  model  this  experiment  using  our  inviscid  code,  we  had  to 
first  calibrate  it  to  achieve  the  same  time  of  arrival.  We  found  that 
increasing  the  initial  conditions  within  the  driver  region  by  5%  was  sufficient 
to  achieve  the  correct  time  of  arrival.  Therefore,  the  following  ICs  were  used 
to  model  the  shock  tube: 

Wi  =[4:J10kg  /  ,0.0  m/s,  398. 1x10^  Pa)  x<0.91m 

W(x,0)  =  j  ^  ^  ;  .  ( 107 ) 

W4  =  (l.208feg/m^,O.Om/s,102.1xlO^Paj  x>0.91m 

To  model  the  left  shock  tube  wall,  we  placed  a  reflective  boundary  at 
X  =  Om .  On  the  left,  we  restricted  the  computational  domain  by  placing  a 
transmissive  boundary  at  x  =  50m .  This  reduced  the  execution  time  required 
but  was  far  enough  from  the  measurement  station  not  to  affect  the  results. 
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The  restricted  spatial  domain,  [Om,  50/n] ,  was  discretized  into  5000  cells. 
Cmax  was  set  to  0.9  and  we  used  the  van  Leer  A  weight  limiter  function. 

Figure  27  shows  the  overpressure  time  history  obtained  via  our  code. 
Unlike  the  Army  Research  Laboratory  code  results,  the  overpressure  history 
curve  is  very  smooth  near  the  shock  front.  This  illustrates  the  impact  our 
TVD  weight  hmiter  functions. 


Figure  27:  Overpressure  Time  History  at  31.44m 
The  shock  arrived  at  our  computational  measurement  station,  located 
at  a;  =  31.44  m,  after  66.08/ns.  The  code  predicts  a  peak  overpressure  at  the 
shock  front  of  61.8fePa.  This  is  within  6.8%  of  the  experimental  results.  This 
is  quite  reasonable  and  suggests  our  code  models  shock  propagation  correctly. 
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Chapter  5:  Shock  Code  Implementation  in  Two  Dimensions 

A.  Initial  Boundary  Value  Problem 

In  this  chapter  we  consider  and  solve  the  two-dimensional  IBVP; 


IBVP  = 


Ut+F(U),+G(U)y^0 

U(x,y,0)  =  U^(x,y) 


(108) 


where  U^{x,y)  is  a  piecewise  function  that  defines  initial  conditions  over  the 
rectangular  spatial  domain,  [0,L^]x  [o,L^J,  at  t  =  Os.  In  two  dimensions  we 

describe  physical  conditions  along  all  four  boundaries  again  using  reflective 
and  transmissive  boundary  conditions. 

B.  Discretization  of  Domain  with  Two  Spatial  Dimensions 

To  discretize  the  domain,  [0,L^]x[0,L^]x[0,r],  we  simply  extend  our 

one-dimensional  mesh  in  the  y  direction.  With  an  additional  spatial 
direction  we  now  have  three  mesh  indices:  i ,  j  and  n  where  the  cell- 
averaged  value  in  cell  (i,j)  at  is 


U^j^U(Xi,yjr) 


(109) 


The  spatial  domain,  [0,L3|.]x  [0,L^J,  is  subdivided  into  I  xj  cells  of 
uniform  dimension,  AxxAy,  where 


Ax  = 


(110) 


and 
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(Ill) 


The  center  of  cell  (i,  j)  is  located  at: 


(112) 


The  physical  extent  of  cell  (i,j)  is  defined  by  its  faces.  The  location  of  each 
face  is  defined  by  its  corresponding  plane  as  follows: 


and 


face  i-1/2:  Xi_y2  = 

(113) 

face  i +  1/2: 

(114) 

f  -I/O  L  (j-1) 

face  j  - 1/2 :  yj_y2  =  - , 

d 

(115) 

face  j  + 1/2 :  y^+ j/2  =  . 

d 

(116) 

The  temporal  domain  is  subdivided  adaptively  into  N  time  steps  of 
non-uniform  length,  At .  Here  we  concern  ourselves  with  limiting  wave 
propagation  in  both  spatial  directions  to  maintain  the  Courant-Friedrichs- 
Lewy  stability  requirement.  In  two  dimensions  therefore,  we  calculate  the 
71*1'  time  step  in  the  following  manner  (Toro,  1997): 

r  ^ 


=  C  min 


A:!c 


Ay 


WL^’W), 


(117) 


max  / 


where  the  maximum  wave  speed  estimates  in  each  spatial  direction  are 
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Wlmx  =°iax(|u,j|+qj)  (118) 

and 

fe'L,  ( 119) 

^yJ 

The  Courant  coefficient,  C,  is  introduced  in  Equation  ( 117  )  to  ensure 
stability  as  described  in  Section  E  of  Chapter  4. 

C.  Dimensional  Splitting 

We  solve  our  two-dimensional  IBVP  approximately  using  a  numerical 
technique  called  dimensional  splitting.  We  choose  this  because  its 
implementation  is  straightforward,  building  directly  upon  our  one¬ 
dimensional  code.  In  dimensional  splitting,  we  split  the  original  IBVP  into 
two  one-dimensional  IBVPs 


ZIBVP  = 


'Ut+F(U)^=0 

U(x,y,t) 


(120) 


and 


yiBVP  = 


'Ut+G(U)y-0 

U(x,y,t) 


(121) 


In  implementing  dimensional  splitting,  our  goal  is  to  develop  a  second- 
order  accurate  solution  procedure  that  remains  computationally  efficient.  To 
do  so,  we  use  the  Beam  and  Warming  variant  (1976): 

^n+2  =  ( 122  ) 
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is  the  X  -flow  time  step  operator  that  operates  on  the  initial  data,  U-^j , 

to  obtain  an  interim  solution,  ,  by  solving  J,  one -dimensional 

problems  in  the  x  -direction  .  Similarly,  we  define  as  the  y  -flow  time 
step  operator. 

To  illustrate  how  dimensional  splitting  is  implemented  via  Beam  and 
Warming,  we  consider  the  following  scenario.  Given  Ui  j ,  Equation  (  50 )  is 

solved  over  the  interval  to  obtain  the  second  order-accurate 

solution,  .  Before  we  can  begin  our  one-dimensional  sweeps  we  must 

first  set  our  boundary  conditions  and  calculate  a  stable  time  step.  Note  from 
Equation  ( 122  )  that  At  remains  constant  over  the  entire  time  interval.  This 
is  necessary  so  that  first-order  errors  will  cancel,  leaving  only  second-order 
errors. 

Once  we  have  determined  At ,  we  begin  the  first  of  two  full  time  steps 
with  .  During  this  step,  X^^  operates  on  U-^j  by  solving  J  one¬ 
dimensional  problems  in  the  x  direction.  Once  the  interim  solution,  , 

is  obtained,  operates  on  it  solving  I  one-dimensional  problems  in  the 
y  direction  to  obtain  .  At  this  point  we  have  completed  one  full  time 

step.  As  a  precautionary  measure,  we  verify  that  the  Courant-Friedrichs- 
Lewy  stability  requirement  has  been  met.  If  our  choice  of  At  was 
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inappropriate,  we  update  C ,  throw  away  our  intermediate  calculations  and 
re-execute  the  first  time  sweep. 

Once  a  stable  solution  is  obtained  at  ,  we  begin  the  second  time 
step  over  the  interval  .  Because  the  method  employed  during  the 

first  time  step  results  in  an  interim  solution  that  is  biased  in  the  y  direction, 
we  now  execute  the  second  time  step  in  reverse  order  to  obtain  the  solution, 

i.e.  operates  on  U-^j^  followed  by  X^‘on  to  obtain  Given 

the  symmetry  of  Equation  ( 122  )  and  a  constant  At ,  this  results  in  a 
cancellation  of  the  first-order  error  terms  giving  us  the  second-order  accuracy 
we  desire.  Before  we  move  forward  in  time,  we  again  verify  the  stability  of 

the  solution.  Once  we  are  satisfied  that  Ui^J^  is  the  proper  solution,  we 

update  our  simulation  clock  by  2  At . 

The  dimensional  splitting  technique  is  repeated  as  many  times  as 
required  to  obtain  the  solution  at  the  desired  time  as  shown  in  Figure  28.  It 
is  important  to  note  that  the  dimensional  splitting  technique  used  here  is 
inherently  parallel.  Calculations  performed  during  each  time  step  operator 
may  be  executed  independently.  With  future  parallel  implementation  in 
mind,  we  wrote  the  code  using  Fortran  90  and  Fortran  95  language 
extensions.  Subroutines  and  functions  were  written  to  be  pure  and 
independent.  To  make  the  code  parallel,  all  that  remains  is  to  add  the 
required  High  Performace  Fortran  language  extensions. 
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While  this  method  is  simple,  the  two-dimensional  nature  of  our  time 
step  operators  complicates  the  solution  procedure.  In  the  next  section,  we 
address  these  special  considerations. 


Set  t  =  Q 

Set  initial  conditions 
Do  until  t>T 


^old 

Set 


=  U 
=  W 


Set  boundary  conditions  in  each  phantom  cell  (in  order  to) 
Calculate  maximum  wave  speeds  in  x  and  y  directions 


Do  time  step  calculations  until  a  stable  second-order  solution  is  obtained 
Calculate  At 

Do  X  direction  one-dimensional  problems  independently  for  j  =  \  to  J 

For  all  cell  interfaces,  t  =  -1  to  I  +  1,  solve  the  Riemann  problems 
For  all  cell  interfaces,  i  =  0  to  I ,  calculate  the  limited  WAF 
For  all  cells,  i  =  1  to  I ,  update  U  and  W 

End  Do 


Set  boundary  conditions  in  y  phantom  cells 

Do  y  direction  one-dimensional  problems  independently  for  i  =  1  to  7 

For  all  cell  interfaces,  j  =  -l  to  J  +  \,  solve  the  Riemann  problems 
For  all  cell  interfaces,  j  =  0  to  J,  calculate  the  limited  WAF 
For  all  cells,  7=1  to  J ,  update  tJ  and  W 

End  Do 

Figure  28:  Pseudocode  for  2D  Code 
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If  fmax  <  1  first  time  step  calculation  is  stable 
Set  boundary  conditions  in  y  phantom  cells 

Do  y  direction  one-dimensional  problems  independently  for  i  =  1  to  / 

For  all  cell  interfaces,  j  =  -l  to  J  + 1 ,  solve  the  Riemann  problems 
For  all  cell  interfaces,  j  =  0  to  J ,  calculate  the  limited  WAF 
For  all  cells,  j  =  1  to  J ,  update  U  and  W 

End  Do 

Set  boundary  conditions  in  x  phantom  cells 

Do  X  direction  one-dimensional  problems  independently  for  j  =  1  to  J 

For  all  cell  interfaces,  i  =  — 1  to  /  +  !,  solve  the  Riemann  problems 
For  all  cell  interfaces,  i  =  0  to  I,  calculate  the  limited  WAF 
For  all  cells,  i  =  1  to  I ,  update  U  and  W 

End  Do 

If  i^niax  <  1 second  time  step  calculation  is  stable 


Exit  second-order  time  step  loop 
Else 


Update  Courant  coefficient 
Set  U  =  Uqi^ 

Set  W  =  W,ia 

Start  second-order  time  step  over 
End  If 
Else 


Update  Courant  coefficient 
Set  U  =  Uqi^ 


Set  W  =  W,ia 

Start  second-order  time  step  over 


End  If 


End  Do 

Update  simulation  time:  t  =  t  +  2At 


End  Do 


Figure  28  continued 
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D.  Special  Considerations 

Recall  from  Equation  (  8 ),  that  in  two  dimensions  U  is  a  four- 
component  vector:  p  (mass),  pu  momentum),  pv  (y  momentum)  and  E 

(total  energy).  Therefore,  special  care  must  be  taken  in  solving  our  one¬ 
dimensional  sweeps.  To  remain  directionally  independent,  we  define  W  and 
U  in  terms  of  normal  and  tangential  components,  i.e. 

W  =  {p,V^,Vt,p)  (123) 

and 

U  =  {p,P^,Pj,,E)  (124) 

where  ,  Vy ,  ,  and  Py  represent  the  normal  and  tangential 

components  of  velocity  and  momentum  respectively.  So,  for  the  general  one¬ 
dimensional  problem,  we  solve  a  system  of  four  PDEs. 

Recall  that,  given  m  PDEs,  we  have  m  waves  in  the  Riemann  problem 
structure.  Therefore,  in  two  dimensions,  the  general  Riemann  problem  wave 
structure  in  the  one-dimensional  WAF  consists  of  four  waves:  wave  1  (left 
outer  wave),  wave  2a  (contact  surface),  wave  2b  (tangential  velocity  shear 
wave)  and  wave  3  (right  outer  wave)  as  shown  in  Figure  29.  We  denote  the 
inner  waves  as  wave  2a  and  wave  2b  to  reinforce  the  fact  that  they  travel  at 
the  same  speed  and  are  therefore  indistinguishable  in  the  x-t  plane. 
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Figure  29:  WAF  Wave  Structure  in  Two-Dimensions 
When  we  solve  the  Riemann  problem  via  Godunov’s  method,  the 
solution  procedure  for  the  tangential  components  is  trivial.  When  we  apply 
limited  WAF  however,  the  solution  procedure  becomes  more  complicated. 
Recall  that  in  limited  WAF,  each  wave,  Wf,  ,  is  amplified  by  its  corresponding 

wave  amplifier,  ,  which  is  dependent  upon  .  T5Tpically,  r2a  *  r2{, .  This 

results  in  a  dispersed  Riemann  problem  wave  structure  that  complicates  the 
solution  procedure.  In  Figure  30,  we  see  that  the  application  of  limited  WAF 
results  in  four  waves  and  five  regions. 
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Figure  30:  Limited  WAF  Wave  Structure  in  Two-Dimensions 
We  calculate  the  limited  weighted  average  flux,  ,  as  a  two  step 

process  where  the  normal  and  tangential  components  are  calculated 
separately.  To  illustrate  how  this  is  done,  we  consider  a  one-dimensional 
calculation  in  the  x  direction.  As  described  in  Chapter  4,  we  calculate 

pLWAF  weighted  average  state,  .  To  calculate  , 

^N^+y2  simply  ignore  the  tangential  velocity  shear 

wave,  W2I) ,  and  solve  the  problem  as  previously  described.  The  flux  of  p  and 
pV^  is  then  calculated  directly. 

Notice  that  p  Vj<  and  E ,  and  their  associated  fluxes,  are  dependent 
upon  the  tangential  velocity  component,  V'p .  We  calculate  V^Y^2  solving 
a  hypothetical  Riemann  problem  separately  that  consists  of  a  single 
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wave, wave  2b  ,  and  two  constant  state  regions,  and  V7-  4 .  The  tangential 
velocity  is  calculated  by  (Toro,  1997): 

^tJ^2  •  ( 125  ) 

The  limited  weights,  W5  and  Wg ,  are  defined  as  (Toro,  1997): 

W5=|(l  +  ^>2b)  (126) 

and 


We  =-|(l-<^2b) 


(127) 


where  -  ^2b  ^2b  •  Unlike  our  typical  TVD  weight  limiter  functions 
(A  >  ^2a  >  ^3)’  ^2b  dependent  upon  the  tangential  velocity  gradient  and  not 
the  density  gradient,  i.e. 


^^Ti+lf2-a,2b 

^i+l/2,2b  -~rT? - 

^^Tuy2,2b 


(128) 


We  now  have  a  solution  procedure  for  every  component  of  and 


pLWAF 


No  other  special  considerations  are  required. 


E.  Sod  Test  in  Two-Dimensions 

Before  we  look  at  more  complex  problems,  we  solved  the  Sod  problem 
in  two-dimensions  to  test  our  code.  The  two-dimensional  Sod  test  was 
performed  over  a  square  computational  spatial  domain,  [0,1]  x  [0,1] .  To 

maintain  the  one-dimensional  nature  of  the  problem  we  apphed  reflective 
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boundaries  at  =  0  and  y  =  1.  Each  parameter  was  set  to  the  same  values 
as  before.  We  applied  the  van  Leer  A  TVD  weight  limiter  function  to  achieve 
monotonic-like  results  near  discontinuities.  The  following  initial  conditions 
were  used: 

ICs: 

3/ 0)  =  r  ^  " (l.0%//n^  Q.Om/s,  O.Om/s, l.OPaf  x  <  0.5/n  ( 129 ) 

[W4  =(o.l25/j^/m^,O.Om/s,O.Oni/s, 0.10Pa|^  x>0.5m 

Figure  31  shows  our  two-dimensional  results.  These  are  in  agreement  with 
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Figure  31:  2D  Sod  Density  Profile  at  t=0.25s 
those  obtained  in  one  dimension.  In  Sod’s  two-dimensional  problem,  flows 
occur  in  only  one  direction.  To  ensure  the  code  models  two  component  flows 
correctly,  we  solved  a  series  of  cyhndrical  explosion  problems  as  presented  by 
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Toro  (1997).  In  each  case,  we  placed  a  high-pressure  cylindrical  region  at  the 
center  of  a  square  computational  domain.  The  solution  of  each  was  verified 
in  the  following  manner: 

-  primitive  variable  profiles  were  qualitatively  compared  against 
Toro’s  results  at  planes  G  =  n  as  measured  from 

the  axis  to  the  leading  edge  of  the  shock  front, 
the  cylindrical  symmetry  of  the  solution  was  verified. 


87 


Chapter  6:  High  Pressure  Cylinder  Problem 

In  Chapter  5  we  extended  our  one-dimensional  shock  code  into  two 
dimensions  using  dimensional  splitting  but  only  gave  a  cursory  look  at  the 
numerical  capabilities  of  the  code.  Here  we  further  investigate  the  numerical 
properties  and  capabilities  of  the  code  to  solve  complex  problems  of  interest  to 
DoD;  namely,  cylindrical  shock  wave  propagation  and  reflection  over  ideal 
surfaces. 

A.  Problem 

To  illustrate  the  numerical  capabihties  of  our  two-dimensional  shock 
code,  we  simulate  the  complicated  shock  phenomena  that  occur  from  the 
detonation  of  a  cylindrically  symmetric  explosive  4.0m  above  an  ideal  surface. 
In  our  model,  we  completely  ignore  the  detonation  process  itself  and  concern 
ourselves  with  its  mechanical  effects  only;  the  outward  propagation  of  the 
resulting  cylindrical  shock  wave.  To  simulate  the  near  instantaneous 
deposition  of  energy  that  occurs  and  the  resulting  high  pressures 
immediately  after  detonation,  we  introduce  a  hypothetical  cylinder  normal  to 
the  x-z  plane  centered  at  (x,z)  =  (Om, 4m)  with  r  =  0.25m .  The  ambient  air 
within  the  cylinder  is  heated  until  we  reach  three  times  the  ambient 
pressure,  i.e.  p^yi  =  3  pg  •  This  gives  us  an  infinitely  long  high-pressure 

cylinder  parallel  to  the  ideal  surface  that: 
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is  the  physical  analogue  of  an  annular  shock  tube,  i.e.  removal  of 
the  cylindrical  diaphragm  begins  the  simulation. 

-  produces  a  cylindrical  shock  wave  shortly  after  removal  of  the 
diaphragm  similar  to  that  produced  by  the  explosive. 

-  is  simple  to  model  numerically. 

B.  Shock  Wave  Propagation  and  Reflection  over  an  Ideal  Surface 

Before  solving  the  problem  numerically,  we  discuss  the  behavior  of  the 
shock  wave  as  it  propagates  outward  and  reflects  over  an  ideal  surface.  The 
qualitative  discussion  below  provides  the  necessary  background  to 
understand  and  interpret  the  numerical  results  presented  in  Section  C.  For  a 
more  detailed  discussion  of  shock  reflection  phenomena  we  suggest  the 
following  books:  Gabi  Ben-Dor’s  Shock  Wave  Reflection  Phenomena.  Samuel 
Glasstone  and  Philip  Dolan’s  The  Effects  of  Nuclear  Weapons  and  Gilbert 
Kinney’s  Explosive  Shocks  in  Air. 

Shock  Wave  Propagation 

We  begin  our  discussion  just  after  the  detonation  of  the  explosive  at 
f  =  Os  as  modeled  by  our  high-pressure  cylinder  approximation.  The 
simulation  begins  with  the  instantaneous  removal  of  the  cylindrical 
diaphragm.  With  the  removal  of  our  fictitious  barrier,  the  heated  air  is  no 
longer  constrained  and  begins  to  accelerate  rapidly  outward  down  the 
pressure  gradient.  When  the  particle  velocity  of  the  heated  gas  exceeds  the 
ambient  speed  of  sound,  a  cylindrical  shock  front  is  formed  that  propagates 
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radially  outward.  As  the  shock  front  continues  its  outward  expansion, 
cylindrical  divergence  of  the  flow  results  in  a  decreasing  peak  overpressure  at 
the  front,  i.e.  oc  1/r . 

Directly  behind  the  shock  front,  there  is  an  outwardly  expanding 
cylindrical  contact  discontinuity.  As  in  the  one-dimensional  case,  this  is  due 
to  the  difference  in  internal  energy  between  the  shock-heated  air  ahead  of  the 
contact  discontinuity  with  the  cooler  expanded  air  behind  it.  In  our 
investigation  of  shock  propagation  and  reflection,  the  contact  discontinuity  is  . 
of  little  interest.  We  simply  mention  it  here  for  completeness. 

At  the  same  time  that  the  shock  front  and  contact  discontinuity  are 
formed,  a  rarefaction  wave  begins  to  travel  up  the  pressure  gradient  towards 
the  axis  of  the  high-pressure  cyhnder.  At  the  instant  the  wave  reaches  the 
axis,  it  instantaneously  reflects  the  cylindrical  wave.  The  resulting 
rarefaction  wave  accelerates  as  it  travels  outward  through  the  previously 
shock-heated  air.  Eventually  the  rarefaction  head  catches  up  to  and  interacts 
with  the  shock  front  resulting  in  the  characteristic  coupled  wave  structure 
shown  in  Figure  32. 
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Figure  32:  Typical  Overpressure  Profile 
Notice  that  the  decrease  in  pressure  behind  the  shock  front  is  due  entirely  to 
the  interaction  of  the  reflected  rarefaction  wave.  In  addition,  the  minimum 
overpressure  within  the  rarefaction  becomes  negative  at  some  point  in  time. 

The  cylindrical  shock  wave  above  continues  to  propagate  as  described 
until  ambient  conditions  are  reached.  Near  the  surface  however,  the  flow  is 
far  more  interesting  and  complicated.  When  our  initial  shock  wave  reaches 
the  surface,  various  reflection  phenomena  occur.  We  will  observe  the 
following  shock  reflection  phenomena:  normal  reflection,  regular  reflection, 
transition  from  regular  to  Mach  reflection,  and  Mach  reflection. 

Shock  Wave  Reflection 

Up  to  now  we  have  concerned  ourselves  with  the  formation  and 
propagation  of  the  shock  wave  over  the  temporal  interval  as  shown  in 

Figure  33.  At  1  =  1^,  the  incident  shock  wave,  I ,  is  traveling  perpendicular  to 
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Figure  33:  Typical  Shock  Reflection  Phenomena 


the  surface  and  normal  reflection  occurs  instantaneously  at  the  point  x  =  Om. 
When  the  flow  is  normal  to  the  surface  meaning  no  tangential  components 
are  present,  the  resulting  peak  overpressure,  ,  is  at  its  maximum 

possible  value,  Ap^ .  The  reflected  peak  overpressure  at  normal  reflection, 

Ap^ ,  is  a  least  twice  the  incident  overpressure,  Apf^^  (Glasstone  and  Dolan, 
1977); 


Ap^=2Apr" 


7po+4Apr^ 

7po+Apr* 


where  po  is  the  ambient  pressure  before  detonation. 


(130) 


As  the  incident  shock  wave  continues  to  propagate,  e.g.  <t<tc,  its 

direction  of  travel  becomes  oblique  with  respect  to  the  surface.  Immediately 
after  t  ,  a  reflected  shock  wave  (i?)  is  formed  marking  the  end  of  normal 

reflection.  The  resulting  two-shock  configuration,  diagramed  in  Figure  34,  is 
called  regular  reflection. 
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Figure  34:  Regular  Reflection 

With  reference  to  Figure  34,  P  is  the  reflection  point  where  wave 

I  and  R  intersect  at  the  surface.  The  peak  overpressure,  ,  occurs  at 

this  point  due  to  the  stagnation  of  the  normal  components  of  the  flow  behind 
I  as  reflection  occurs.  As  /  continues  to  propagate  over  the  surface,  the 
angle  of  incidence,  0j ,  increases  resulting  in  a  corresponding  decrease  in  the 

normal  flow  component.  In  addition,  Apf^^  decreases  as  described  earlier 

due  to  cylindrical  divergence.  Therefore  decreases  and  is  always  less 

than  Ap^. 

Regular  reflection  continues  until  a  critical  angle,  ,  is  achieved.  In 
Figure  33,  this  occurs  at  t  •  At  this  point,  we  begin  the  transition  from  a 

two-shock  regular  reflection  to  a  three-shock  configuration  called  Mach 
reflection.  This  occurs  as  R  begins  to  coalesce  into  I  forming  a  single  shock 
front  called  the  Mach  stem  (M).  This  occurs  as  the  shock  moves  through  the 

f 

transition  region,  during  the  interval  <t  <1^.  . 
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Mach  reflection  begins  when  the  Mach  stem  is  formed  as  depicted  in 
Figure  35.  The  peak  overpressure  at  the  Mach  stem  is  greater  than  that 

observed  for  either  /or  i?,  i.e.  >  Ap“^* .  As  we  move  up  M , 

the  overpressure  decreases.  The  intersection  of  M ,  R  and  I  is  called  the 
triple  point  (T).  Notice  in  Figure  33  that  for  t>tc  ,  the  height  of  T  increases 
over  time  as  M  continues  to  grow. 


Now  that  we  have  an  understanding  of  shock  propagation  and 
reflection,  we  present  and  discuss  the  numerical  results  obtained  using  our 
two-dimensional  shock  code  in  the  next  section.  If  our  two-dimensional  code 
models  the  physics  of  inviscid  flow  in  two-dimensions  correctly  and  our 
implementation  is  sound,  we  should  observe  the  shock  phenomena  described 
above. 

C.  Numerical  Solution 

To  model  the  complicated  shock  phenomena  that  results  from  the 
detonation  of  our  cylindrical  explosive  we  solved  the  two-dimensional  IBVP 
given  the  ICs  in  Equation  ( 131 )  and  BCs  in  Table  5.  Recall  the  high- 
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(131) 


=(l.2045fe^/m^0.0m/s,0.0m/s,304.05xl0^Paf  r<0.25m 
^4  =  (1.2045 %/m^,0.0m/s,0.07n/s,101.35xl0^Paf  r> 0.25m 


X  BCs 

2  BCs 

Reflective  along  x  =  Om . 

Reflective  along  z  =  Om . 

Transmissive  along  x  =  20m  . 

Transmissive  along  z  =  12m . 

Table  5:  High-pressure  Cylinder  Boundary  Conditions 
pressure  cylinder  models  the  energy  deposition  immediately  after  detonation 
and  is  centered  at  {x,z)  =  (Om,  4m) .  To  reduce  the  execution  time  of  the  code, 

we  took  advantage  of  the  symmetry  of  the  problem  and  inserted  a  reflective 
boundary  along  the  x  =  Om  plane.  Our  remaining  reflective  boundary 
simulated  the  ideal  surface,  i.e.  the  ground.  To  restrict  our  computational 
domain  to  a  finite  size,  we  employed  transmissive  boundaries  at  x  =  20m  and 
z  =  12m .  Each  was  positioned  at  a  sufficient  distance  to  allow  observation  of 
interesting  flow  features  while  delaying  the  arrival  of  the  inevitable  false 
reflected  waves  until  after  the  passage  of  the  shock  through  the  region  of 
interest.  These  boundary  conditions  therefore  define  the  following  spatial 
domain:  [Om,  20m]  x  [Om,  12m] . 

Initially,  we  subdivided  our  spatial  domain  into  600x360  square 
computational  cells.  The  square  cell  geometry  was  used,  instead  of  the  more 
typical  and  computationally  efficient  rectangular  cell,  to  limit  the  numerical 
artifacts  that  occur  in  approximating  our  cylinder  on  a  Cartesian  mesh.  In 
addition,  we  used  a  weighted  average  of  the  initial  conditions  in  diaphragm 
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cells.  As  we  shall  see  in  our  results,  this  significantly  decreased  the  stair¬ 
step-like  numerical  artifacts  that  occur  at  early  times  resulting  in  a  nice 
cylindrical  wave  structure. 


Lastly,  we  set  the  execution  parameters  as  shown  in  Table  6. 


^0 

n 

'-'max 

Qs 

Tolerances 

TVD  Weight  Limiter 

0.2 

0.9 

2 

1.0x10"^ 

van  Leer  A 

Table  6:  High-pressure  Cylinder  Execution  Parameters 
Notice  the  Courant  coefficient  is  set  as  described  in  Chapter  4  to  maintain 
stability  at  early  times  and  ^>vla  was  chosen  to  ensure  sharp  resolution  of  the 
shocks  of  interest. 

We  ran  the  simulation  over  the  time  interval  [Oms,  60ms] .  To  ensure 
we  captured  regular  reflection,  transition  from  regular  reflection  to  Mach 
reflection,  and  Mach  reflection  as  described  in  Section  B,  we  recorded  the 
location  and  value  of  the  maximum  overpressure  every  other  time  step  as 
shown  in  Figure  36-  38.  In  addition,  the  overpressure  solution  was  saved 
every  0.4ms .  This  data  was  then  used  to  visuahze  the  structure  of  the  shock 
waves  at  times  of  interest  using  Tecplot  (ver  7.5),  a  high-end  scientific 
visualization  software  package. 

During  our  discussion  of  the  numerical  results  we  will  follow  the 
maximum  overpressure  curves  shown  in  Figure  36,  37  and  38.  These  curves 
serve  to  guide  our  discussion  and  show  the  progress  of  the  simulation  as  the 
cylindrical  shock  wave  forms,  propagates  outwards  and  reflects  over  the 
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surface.  Periodically,  we  will  look  at  the  overpressure  solution  obtained  at 
various  times  as  weU. 

In  Figure  36,  it  is  immediately  obvious  from  our  time  plot  that  two 

time  intervals  exist:  a  free-air  time  interval  where  the  shock  has  not  yet 

reached  the  surface  and  an  interval  over  which  reflection  occurs.  During  the 

free-air  interval,  [0ms,  9.84ms) ,  the  maximum  overpressure  occurs  along  the 

shock  front  of  the  cylindrical  shock  wave  as  it  propagates  outwards.  We 

found  that  occurs  at  the  shock  front  and  decreases  like  due  to 

/r 

cylindrical  divergence. 


0.0  5.0  10.0  15.0  20.0  25.0  30.0  35.0  40.0  45.0  50.0  55.0  60.0 

t  (ms) 


Figure  36:  Maximum  Overpressure  vs  Time  [0ms,  60ms] 
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In  Figure  39,  we  see  that  the  high-pressure  cylinder  approximation  has 
produced  a  nearly  cylindrical  shock  wave  as  desired.  The  stair-case-like  wave 
patterns  present  along  the  shock  front  are  numerical  artifacts  of  our 
weighted  cylinder  approximation  and  cannot  be  helped.  A  finer  mesh  will 
limit  this  effect  considerably  but  at  the  significant  cost  of  additional 
computational  effort. 

At  t  s  9.847ns ,  we  observed  a  sharp  rise  in  the  maximum  overpressure 
fi:om  \A.2QkPa  to  25.35kPa  as  shown  in  Figure  36.  This  indicates  that  the 
cylindrical  shock  wave  has  reached  the  surface  and  reflected  producing  the 
peak  overpressure  observed  at  the  reflection  point  on  the  surface.  Note  that 
this  rise  in  overpressure  is  inconsistent  with  Equation  ( 130 )  because  we  did 
not  capture  the  exact  moment  at  which  normal  reflection  occurred. 

Figure  37  illustrates  the  maximum  overpressure  versus  time  over  the 
time  interval  [9.847ns,  GO/ns] .  During  this  time  interval  reflection  occurs  and 

the  maximum  overpressure  remains  located  at  the  surface.  At  early  times 
where  0j  <  0"''* ,  regular  reflection  occurs  as  shown  in  Figure  40.  Here  we 
see  an  excellent  example  of  regular  reflection.  Both  the  reflected  and 
incident  waves  are  well  resolved.  In  addition,  we  note  that  the  peak  reflected 
overpressure  occurs  due  to  the  stagnation  of  the  incident  flow  at  the 
reflection  point. 
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Figure  37:  Maximum  Overpressure  vs  Time  [9.84ms,  60ms] 
Notice  that  if  we  extend  the  left  and  right  portions  of  the  curve  as 
depicted  by  the  dashed  lines  in  Figure  37  we  can  clearly  see  a  bump  that 
occurs  roughly  between  10.8  and  13.1ms .  Over  this  time  interval,  the 

behavior  of  changes;  it  decreases  at  a  rate  slower  than  that  observed 

at  earlier  times  until  t  =  13.1ms .  This  change  is  due  to  the  fact  that  the 
reflected  wave  is  catching  up  to  and  stagnating  off  the  incident  wave,  i.e. 
transition  is  occurring  from  regular  to  Mach  reflection.  At  t  =  13.1ms , 

transition  ends  with  the  formation  of  the  Mach  stem  and  decreases  at 

a  faster  rate.  Therefore,  the  reflection  time  interval  may  be  subdivided  into 
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three  distinct  intervals:  [9.84nis  10.8/ns)  -  regular  reflection,  [10.8/ns,  13.1/ns] 
-  transition  and  (13.1/ns,  60/ns]  -  Mach  reflection. 

In  Figure  38,  we  plotted  the  maximum  overpressure  during  reflection 
versus  the  ground  distance  at  which  it  occurred.  Both  are  functions  of  t  and 
give  us  the  time  of  arrival  of  the  shock  on  the  ground.  We  subdivided  the 
reflection  interval  into  regions  just  as  we  did  in  Figure  37.  Here  we  see  that 
the  transition  region  begins  near  x  =  1.68/n  and  ends  around  x  =  3.28/n . 


Figure  38:  Maximum  Overpressure  vs  Ground  Distance 
To  illustrate  and  investigate  the  transition  from  regular  to  Mach 
reflection,  we  refined  the  mesh  by  a  factor  of  two  and  obtained  the 
overpressure  solutions  every  0.4/nsover  the  interval  [10.4/n.s,  13.6/ns]  as 


100 


shown  in  Figxn’e  42  -  49.  As  we  proceed  in  time  through  the  transition  region, 
we  see  that  as  Oj  increases  the  distance  at  the  surface  between  the  reflected 
and  incident  shock  waves  decreases  steadily.  In  addition,  the  width  of  the 
peak  overpressure  region  decreases  but  grows  in  height  suggesting  that  flows 
behind  the  reflected  wave  are  beginning  to  stagnate  as  it  catches  up  to  the 
incident  wave.  This  indicates  that  the  reflected  and  incident  shock  waves  are 
beginning  to  merge  near  the  surface.  By  f  =  12.8ms ,  we  see  the  beginnings  of 
a  Mach  stem  forming  near  the  ground  and  hy  t  =  13.6ms  a  Mach  stem  is 
present. 

In  Figure  51  —  52  we  see  the  well-resolved  three-shock  configuration  of 
Mach  reflection.  As  time  goes  on  and  the  waves  continually  merge  we 
observe  the  growth  of  the  Mach  stem,  the  motion  of  the  triple  point  away 
from  the  surface  and  the  decrease  in  the  peak  overpressure.  In  addition,  we 
notice  that  as  we  move  up  the  Mach  stem  the  stagnation  pressure  decreases 
as  we  expect. 

Just  before  t  =  58ms  we  see  in  Figure  37,  37  and  38  that  the  peak 
overpressure  at  the  Mach  stem  decreases  suddenly.  This  is  a  numerical 
artifact  of  the  right  transmissive  boundary.  We  observe  the  unexpected  dip 
in  peak  overpressure  because  a  false  reflected  shock  wave  has  traveled  back 
to  the  left  and  perturbed  the  solution.  The  dip  indicates  the  leading  edge  of 
that  false  reflected  shock  wave. 
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From  our  qualitative  analysis  we  see  that  the  two-dimensional  shock 
code  is  capable  of  modeling  the  complicated  shock  phenomena  of  interest  to 
DoD. 
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Figure  39:  High-pressure  Cylinder  Overpressure  at  t=4.8ms 
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Peak  Overpressure  =  27.94  kPa 
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Figure  40:  High-pressure  Cylinder  Overpressure  at  t=10.07ns 
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Peak  Overpressure  =  27.12  kPa 
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Figure  41:  High-pressure  Cylinder  Overpressure  at  t=10.4ms 
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Peak  Overpressure  =  26.49  kPa  H  25.87 
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Figure  42:  High-pressure  Cylinder  Overpressure  at  t=10.8/ns 
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Peak  Overpressure  =  25.96  kPa  ■  25.35 

■  24.74 
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Figure  43:  High-pressure  Cylinder  Overpressure  at  t=11.2ms 
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Peak  Overpressure  =  25.58  kPa 
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Figure  44:  High-pressure  Cylinder  Overpressure  at  t=11.6ms 
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Peak  Overpressure  =  25.21  kPa 
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Figure  45:  High-pressure  Cylinder  Overpressure  at  t=12.0ms 
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Figure  46:  High-pressure  Cylinder  Overpressure  at  t=12.4/ns 
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Peak  Overpressure  =  24.31  kPa 
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Figure  47:  High-pressure  Cylinder  Overpressure  at  t=12.8ms 
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Peak  Overpressure  =  23.91  kPa 
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Figure  48:  High-pressure  Cylinder  Overpressure  at  t=13.2/ns 
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Peak  Overpressure  =  23.38  kPa 
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Figure  49:  High-pressure  Cylinder  Overpressure  at  t=13.6ms 
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Peak  Overpressure  =  13.04  kPa 
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Figure  50:  High-pressure  Cylinder  Overpressure  at  t=27.0/ns 
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Peak  Overpressure  =  11 .75  kPa  ■  11.13 
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Figure  51:  High-pressure  Cylinder  Overpressure  at  t=31.2ms 
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Peak  Overpressure  s  10.51  kPa 
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Figure  52:  High-pressure  Cylinder  Overpressure  at  t=36,0ms 
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Chapter  7:  Conclusion 

The  primary  objective  of  this  effort  has  been  to  develop  a  two- 
dimensional  hydrodynamic  shock  code  based  upon  E.  F.  Toro’s  weighted 
average  flux  (WAF)  method  for  the  investigation  of  air  blast  phenomenology. 
In  support  of  that  effort,  we  initially  developed,  tested  and  validated  a  one¬ 
dimensional  shock  code  based  upon  Toro’s  weighted  average  flux  (WAF) 
method.  During  development  and  testing  of  the  one-dimensional  shock  code, 
we  observed  the  following: 

-  Godunov’s  method  accurately  predicted  the  location  of  all  discontinuities 
but  the  solution  was  smeared  due  to  numerical  viscosity 

-  application  of  WAF  resulted  in 

•  decreased  numerical  viscosity 

•  oscillations  near  discontinuities  typical  of  second-order  accurate 
methods 

-  application  of  TVD  weight  hmiter  functions  adequately  removed 
oscillations  while  maintaining  behavior  indicative  of  second-order 
methods 

-  among  the  weight  limiter  functions  tested, 

•  Super  A  resolved  contact  discontinuities  best 

•  van  Leer  A  weight  limiter  function  resulted  in  smooth  regions  behind 
shocks. 
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Before  moving  into  two  dimensions,  we  validated  the  one-dimensional  shock 
code  against  experimental  results  obtained  at  Army  Research  Laboratory’s 
57cm  shock  tube  facility.  Here,  we  found  the  one-dimensional  shock  code 
underestimated  the  peak  overpressure  at  the  shock  front  by  6.8%. 

Once  the  one-dimensional  shock  code  was  complete,  we  extended  it 
into  two  spatial  dimensions  via  Warming  and  Beam’s  variant  of  dimensional 
spHtting.  To  verify  proper  implementation  of  this  technique,  we  solved  both 
Toro’s  cylindrical  explosion  problem  and  Sod’s  shock  tube  problem  in  two 
dimensions.  Having  shown  the  two-dimensional  code  worked  properly,  we 
modeled  the  detonation  of  a  cylindrically  symmetric  explosive  over  an  ideal 
surface  to  illustrate  the  capabilities  of  the  code.  We  found  the  numerical 
solution  modeled  the  physics  of  shock  propagation  and  reflection  well  and 
that  the  transition  from  regular  to  Mach  reflection  occurred  as  expected.  As 
a  result  of  this  effort,  Air  Force  Institute  of  Technology  now  has  a  two- 
dimensional  hydrodynamic  shock  code  to  further  investigate  air  blast 
phenomenology. 


118 


Appendix  A:  Toro’s  Initial  Guess  Generator 


Rapid  convergence  is  achieved  when  is  sufficiently  close  to  the 

root.  If  is  far  from  the  actual  root,  the  computationally  efficiency  of  ERS 
is  severely  degraded  since  much  of  the  CPU  time  solving  Equation  ( 18  ),  and 
therefore  computational  effort,  is  attributed  to  finding  the  root  of  equation 
(  33 ).  Worse,  an  extremely  poor  choice  can  result  in  a  physically  meaningless 
solution,  p,  <  0.  Clearly,  success  is  dependent  upon  the  numerical  nature  of 

Equation  ( 30 )  and  our  choice  of  pi^^ . 

Toro’s  pressure  function,  f ,  is  monotonic  and  concave  down  because 
/'  >  0  and  /’  <  0 .  This  means /is  relatively  benign.  Therefore,  the  major 
concern  here  is  not  whether  convergence  is  achieved  but  rather,  how  quickly 
it  is  achieved.  Graphical  inspection  of  the  mathematical  properties  of 

Equation  ( 30 )  provides  some  insight  into  how  an  appropriate  p^^  might  be 
chosen  given  the  initial  conditions  and  some  assumptions. 
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ApAA) 


Figure  53:  General  Behavior  of  Toro's  Pressure  Function 
Figure  53  illustrates  the  general  behavior  of  the  pressure  function 
given  initial  conditions,  Wi  and  W4  (Toro,  1997).  Each  curve  in  Figure  53 
represents  the  pressure  function  given  the  same  set  of  initial  conditions  for 
density  and  pressure.  The  difference  between  the  curves  is  due  to  the 
particle  velocity  difference,  where  Aw  =  M4  -  w^  varies. 

With  reference  to  Figure  53,  the  minimum  and  maximum  pressures 
are  defined  to  be: 

Pmin  =min(pi,P4)  (132) 

and 

^max  max(pi,p4).  (133) 
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These  variables  divide  the  pressure  interval  into  three  different  subintervals: 

^2  ■  Pmm  ^  P  —  Pmax 

and 

-^3  •  Pmax  <  P  ^  . 

Note  that  depending  upon  the  value  of  Au ,  the  pressure  function  is 
translated  up  or  down  and  the  root  changes  accordingly. 

For  sufficiently  large  Aw  the  root  lies  within  the  first  interval  such 
thatp,  <  <  Pmax  •  This  indicates  a  two-rarefaction  configuration.  For 

(Au)2  ,  p.  lies  within  the  second  interval  where  <  p.  <  p^ax  •  Then  one 

of  the  outer  waves  is  a  rarefaction  and  the  other  a  shock.  Lastly,  for  small 
Au ,  the  root  lies  within  the  third  interval  where  p,  >  p„,„  >  p.  .  This 

indicates  the  two-shock  configuration. 

Given  the  general  behavior  of  the  pressure  function  and  initial 
conditions,  we  can  predict  the  general  wave  structure  of  the  Riemann 
problem  solution  before  it  is  solved.  Armed  with  this  knowledge,  we  can 

predict  pl®^  sufficiently  close  that  the  converged  solution  is  t5T)ically  obtained 
within  ten  iterations. 

To  determine  where  the  root  lies,  we  begin  with  the  following 
assumptions: 

-  the  root  lies  within  the  second  interval. 
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flow  conditions  are  mild  enough  that  a  linearized  approximation  is 
sufficient. 


With  these  assumptions,  is  obtained  in  the  following  manner  (Toro, 
1991): 


=  max(1.0  X 10  ^Pa,Pp^j.s )  ( 134 ) 

where 


1  1 

Ppvrs  ”■  ~r( Pi  )( Pi  Pa)(^1  ( 135  ) 

2  8 

is  taken  from  Toro’s  PVRS.  Equation  ( 134 )  is  applied  to  ensure  pi®^  >  0 
because  Equation  ( 135  )  may,  incorrectly,  predict  a  negative  pressure.  Note 
that  if  strong  shocks  are  present,  our  linearity  assumption  is  no  longer  valid 
and  PVRS  does  not  perform  well.  Flow  conditions  are  considered  mild 
provided: 


^  Pmax  _ _ _ 

Vs  > -  ( 136 ) 

Pmin 

where  Q^is  a  user-defined  value  typically  set  at  2  (Toro,  1991).  Even  if 

Equation  ( 136 )  is  false,  we  still  calculate  ( 134 )  for  use  in  the  two-shock 
approximation. 

At  this  point,  we  check  where  pl®^  lies  in  relation  to  p^ji^  and  p^gx  •  If 
Pmin  <  -  Pmax  ^^d  flow  Conditions  are  sufficiently  mild,  pl®^  is  used  as 

an  initial  guess  in  Equation  (  33 ).  If  p^®^  <  p^^^ ,  we  assume  the  two- 


122 


rarefaction  configuration  and  is  obtained  using  the  two-rarefaction 
approximation  derived  from  Equation  (  30 )  (Toro,  1997): 


Cl  +C4  --(r-l)(U4  -  Ml) 
+.  ^4 


{r-V)l2r  '  (r-i)/2r 

Pi  Pa  j 


(137) 


Note  that  if  two-rarefactions  are  actually  present  in  the  solution  no  iterations 
are  required  since  this  approximation  is  exact.  If,  however,  >  Pmax^^ 


Qg  <  is  false,  the  two-shock  configuration  is  assumed  and  the  following 


i^min 


two-shock  approximation  is  used  (Toro,  1997): 


Pi 


Pa 


2 


N  1/ 
\/2 


Pi(pi®\r+i)+Pi(r-i) 

2 


\  1/ 

^/2 


P4(p*”^(r+i)+P4(r-i). 


-I-  Ml  -  M4 


(138) 


f  ^  ^ 

1/ 

72 

^  2  ^ 

1/1 

/2 

VPi(p*‘’^(r +i)+Pi(r-i)> 

^  Pa  (P*°^  (r  + 1)  +  Pa  (r  - 1) ; 
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